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ABSTRACT 


A  zonal  Navier-Stokes  model,  developed  by  J.C.  Wu,  is 
installed  and  verified  on  the  NASA  Ames  Cray  X/MP-48 
computer  and  is  used  to  calculate  the  flow  field  about  a 
NACA  0012  airfoil  oscillating  in  pitch.  Surface  pressure 
distributions  and  integrated  lift,  pitching  moment,  and  drag 
coefficients  and  integrated  lift,  pitching  moment,  and  drag 
coefficients  versus  angle  of  attack  are  compared  to  existing 
experimental  data  for  four  cases  and  existing  computational 
data  for  one  case.  These  cases  involve  deep  dynamic  stall 
and  fully  detached  flow  at  and  below  a  freestream  Mach 
number  of  .184.  The  flow  field  about  the  oscillating 
airfoil  is  investigated  through  the  study  of  pressure, 
vorticity,  local  velocity  and  stream  function.  Finally,  the 
effects  of  pitch  rate  on  dynamic  stall  are  investigated. 
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I.  INTRODUCTION 


Dynamic  stall  is  a  phenomenon  that  refers  to  an  airfoil 
delaying  stall  beyond  its  static  stall  angle  due  to  a  rapid 
change  in  angle  of  attack.  Associated  with  this  is  the 
generation  of  a  strong  vortex  that  appears  at  the  leading 
edge  of  the  airfoil  which  expands  as  it  moves  aft,  causing 
large  excursions  in  the  pressure,  pitching  moment  and  lift 
the  airfoil  experiences  (Figure  1,  taken  from  [Ref.  1]). 

Because  the  dynamic  stall  angle  generally  occurs  at  much 
higher  angles  of  attack  than  the  static  stall  angle,  the 
maximum  lift  the  airfoil  generates  can  also  be  much  higher 
than  for  steady  conditions.  Unfortunately  this  is  a 
transient  condition,  as  the  lift  drops  sharply  when  the 
vortex  is  shed  from  the  trailing  edge.  However,  if  the 
mechanisms  that  govern  the  initiation  and  development  of 
this  vortex  and  the  dynamic  stall  phenomenon  can  be 
understood  and  controlled,  it  could  open  the  door  to  much 
more  maneuverable  and  higher  performing  aircraft.  Nature 
provides  a  prime  example  of  what  is  possible  in  the 
dragonfly,  which  uses  vortex  energy  recovery  to  help  achieve 
its  remarkable  maneuverability. 

NASA  currently  has  a  flight  test  program  utilizing 
vortex  generating  leading  edge  slats  to  gain  further  insight 
into  this  fascinating  area.  Parallel  studies  are  underway 


v 


Mi 


% 


g 

I 

5» 


IV 

B 

IV 

If 


vvv-vs  > 


(»)  STATIC  STALL  ANGLE  EXCEEDED 


(b)  FIRST  APPEARANCE  OF  FLOW 
REVERSAL  ON  SURFACE 


(c)  LARGE  EOOIES  APPEAR  IN 
BOUNDARY  LAYER 


(d)  FLOW  REVERSAL  SPREADS  OVER 
MUCH  OF  AIRFOIL  CHORD 


(•)  VORTEX  FORMS  NEAR 
LEADING  EDGE 


((I  LIFT  SLOPE  INCREASES 


(9)  MOMENT  STALL  OCCURS 


(M  LIFT  STALL  BEGINS 


(i)  MAXIMUM  NEGATIVE  MOMENT 


<j)  FULL  STALL 


(k)  BOUNDARY  LAYER  REATTACHES 
FRONT  TO  REAR 


INCIOENCE,  or.  da« 


(I)  RETURN  TO  UNSTALLED  VALUES 


Figure  1.  Dynamic  Stall 
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by  numerous  groups  to  investigate  dynamic  stall  via 
computational  methods.  The  approaches  that  nave  received 
the  most  attention  have  used  Navier-Stokes  modeling,  which 
has  the  flexibility  to  describe  many  flows. 

Navier-Stokes  modeling  has  usually  been  formulated  along 
velocity/pressure  lines.  This  formulation  has  two  major 
problems  [Ref.  2].  The  first  problem  is  caused  by  the  size 
of  the  flow  field,  which  in  mathematical  terms,  is  infinite 


in  extent. 


Boundary  conditions  at  infinity  must  be 


satisfied  and  an  infinite  flow  field  has  to  be  modeled  by  a 
finite  number  of  grid  points.  To  do  this,  the  boundary 
conditions  must  be  previously  known  and  multiple  boundary 
condition  solutions  computed  and  compared  to  the  known 
solution.  Alternatively,  a  coordinate  transformation  can  be 
made  from  the  infinite  flow  field  to  a  finite  one. 
Unfortunately,  this  introduces  a  requirement  for  additional 
computations  of  an  increasingly  complex  nature. 

The  second  problem  is  the  large  number  of  data  points 
required  to  adequately  model  the  flow  field.  Even  for  two 
dimensional  cases  this  can  be  a  significant  impediment. 
Various  grid  spacing  methods  have  been  used  to  concentrate 
data  points  in  the  regions  around  the  airfoil  where  high 
flow  variable  gradients  occur,  but  the  number  of  data  points 
required  is  still  quite  high. 

An  alternative  method  is  to  use  velocity/vorticity 
modeling.  By  taking  advantage  of  certain  features  of 
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velocity  and  vorticity  fields  for  incompressible  viscous 
flow,  certain  conclusions  may  be  reached  [Ref.  3]. 


1.  Vorticity  can  be  neither  created  nor  destroyed  in  the 
interior  of  the  fluid. 

2.  The  total  vorticity  in  the  infinite  unlimited  space 
jointly  occupied  by  the  fluid  and  the  solid  is  always 
zero. 

3.  The  rate  of  convective  transport  of  vorticity  is 

finite.  The  rate  of  diffusive  transport  is 

effectively  finite. 

4.  At  large  distances  from  the  body,  the  velocity  field 
approaches  zero  with  increasing  distance  from  the 
solid. 

5.  At  large  distances  from  the  solid,  the  vorticity  field 
decays  exponentially  with  increasing  distance  from  the 
solid. 

6.  The  velocity  field  is  uniquely  determined  by  the 
vorticity  distributed  in  the  infinite  unlimited  space 
jointly  occupied  by  the  fluid  and  the  solid. 
Alternatively,  the  velocity  field  is  uniquely 
determined  by  the  vorticity  distribution  in  the  fluid 
and  the  velocity  condition  on  the  solid  boundary. 

These  features  have  three  major  salutary  effects.  [Ref.  4] 

1.  The  actual  number  of  grid  points  that  need  to  be 
solved  for  will  generally  be  much  less  than  the  total 
number  of  grid  points  needed  to  define  the  flow  field. 
This  is  because  vorticity  is  generated  solely  at  the 
solid  boundary,  and  this  vorticity  is  diffused  only  a 
short  distance  from  the  solid  before  being  carried 
away  by  convection  and  diffusion.  Therefore,  the  flow 
will  be  inviscid  a  short  distance  ahead  of  the  solid 
and  often  inviscid  at  small  to  moderate  distances 
above  and  below  the  solid. 

2.  By  taking  advantage  of  the  integral  representation 
that  the  vorticity/velocity  equations  can  be  expressed 
in,  the  various  regions  of  the  flow  field  can  be 
treated  separately,  with  separate  grids  and  distinct 
computational  procedures  that  consider  the  length 
scales  and  defining  equations,  with  no  loss  of 
accuracy  and  without  a  need  to  match  the  solutions  of 
the  various  zones.  Ultimately,  an  elegant  method  to 
compute  flow  field  solutions  becomes  available  that 


can  accurately  compute  dynamic  stall  conditions 
without  resorting  to  brute  force  velocity/pressure 
finite  difference  computations. 

3.  The  solution  may  be  expressed  as  an  integral  which  can 
be  converted  to  a  Fourier  series  expansion.  [Ref.  5] 
This  allows  for  very  accurate  computational  solutions 
at  each  grid  point. 

Currently  the  model  is  limited  to  two-dimensional 
incompressible  flows,  but  the  concepts  are  applicable  to 
three-dimensional  and  compressible  flows  as  well. 

A  zonal  Navier-Stokes  model  has  been  developed  by  J.C. 
Wu  and  his  associates  at  the  Georgia  Institute  of  Technology 
and  was  made  available  for  this  study.  An  explicit, 
integro-dif ferential  methodology  procedure  is  used  to  solve 
two-dimensional,  Reynolds  averaged,  incompressible  Navier- 
Stokes  equations. 

The  goals  of  the  present  study  were: 

1.  Install  and  verify  the  code  on  the  Cray  X/MP-48. 

2.  Compare  the  code’s  solutions  with  previous 
experimental  and  theoretical  results. 

3.  Modify  this  code  to  enhance  its  utility. 


II .  MATHEMATICAL  FORMULATION 


A.  GOVERNING  EQUATIONS  DEVELOPMENT 

In  the  study  of  fluid  flows,  normally  certain 
assumptions  will  be  made.  These  assumptions  allow  emphasis 
on  the  specific  features  of  interest,  while  providing  useful 
simplifications  in  the  theoretical  development.  Properly 
chosen,  they  can  also  promote  the  ease  of  numerical 


formulation  and  consequent  solution. 


Accordingly,  the 


assumptions  made  for  the  flow  type  in  this  study  and  their 
respective  purposes  are: 

1.  Two-dimensional.  The  primary  benefit  of  this 

assumption  is  a  very  significant  reduction  in 
computational  requirements.  A  relatively  minor  amount 
of  information  is  lost  for  spanwise  flow  and  tip 
effects,  with  the  salient  aspects  to  be  investigated 
remaining  intact. 

2.  Incompressible.  A  useful  simplification  to  highlight 
the  dynamic  stall  effects  and  allow  for  a  more  readily 
formulated  computational  model. 

3.  Viscous.  A  necessary  element  to  describe  the  solid 
body  and  flow  field  interaction.  This  also  helps  to 
control  the  extent  of  flow  field  computations 
required . 

4.  Unsteady.  Required  to  allow  for  airfoil  pitching. 

5.  Turbulent.  Required  to  describe  the  flow  field  under 
the  conditions  of  interest. 

6.  Reynolds  averaged.  A  formulation  of  the  Navier-Stokes 
equations  for  turbulent  conditions. 


j; 


For  three-dimensional  unsteady  flow  in  the 


direction : 


or 


where : 


Du  _  1  3p  .  .  „2. 


c  pjr  -  pX  -  —  ^  +  yV  u 
Dt  p  dX 


,3u  ,  3u  ,  3u  ,  3u, 

pfcr+Uir-  +  V^-+W^-) 
ot  ax  3y  dZ 


=  MX  -  i  |E  +  +  ifn  + 

p  o  ax2  3y2  Sz2 


^-t^+u^  +  v^  +  w^  =  substantial  derivative  of  u 
Dt  3t  3x  3y  3z 


pX  =  body  forces 


=  pressure  forces 


^2  ,2  .2 

l  ,  3  u  ,  j  u  3  u 


+  - — ^  -t  — =  viscous  term. 


ox  3v  9z 


2 .  Reynolds  Averaging  for  Turbulent  Flows 

The  boundary  layer  equation  of  motion  in  t 
dimensions  is  [Ref.  6]: 


^  =  -  42  +  Uu  £> 


Dt  3x  3y  3y 


letting  the  following  variables  be  defined  as: 


7 


u  =  U  (x,  y)  +  u  '  ( x ,  y  ,  t ) 


v  =  V(x, y)  +  v' (x, y , t) 


p  =  P(x)  +  p'  (x,  y ,  t) 


where : 


U(x,y)  =  mean  x-direction  velocity  value 
over  time 

u'  =  fluctuation  value  of  velocity  in 
x-direction 

u'  «  U 

p  =  pressure 


Then: 


1.  Substituting  u,  v  and  p  into  the  boundary  layer 
equation, 

2.  Taking  the  mean  value  of  each  term, 

3.  Observing  that  the  mean  value  linear  terms  in  a 
fluctuation  component  vanishes, 

4.  Assuming  that  the  derivatives  of  the  mean  value  linear 
term  vanishes,  and 

5.  Neglecting  turbulent  normal  stress  compared  to  the 
shearing  stress, 


gives : 


where : 


,ri  3U  .  „  311 
p(U  3x  V  3y 


dP  ,  3  ,  — ; — 

+  —  -pu'v' 
dx  3v 


-ou'v'  =  the  Reynolds  stress. 
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Generally,  the  neglected  terms  are  much  smaller  than  -pu'v', 
and  it  should  be  noted  that  u'  and  v'  can  reach  values  as 
high  as  .1  U.  However,  for  the  flow  under  consideration, 
these  are  still  valid  approximations. 


Prandtl's  Mixing  Length 


From  the  Reynolds  stress  development,  another 
concept  can  be  presented  [Refs.  7,8].  Boussinesq  introduced 
a  mixing  coefficient  that  related  the  Reynolds  stress  to  the 
derivative  for  U.  This  is  possible  because  for  there  to  be 


a  net  momentum  transfer  from  the  higher  to  lower  momentum 
layers  -pu'v'  >  0  or,  in  other  words,  u'  and  v'  must  be 
positively  correlated.  This  is  done  by  setting  the  Reynolds 
stress  equal  to  the  derivative  of  U  or: 


dU 

dy 


where: 

vT  =  eddy  viscosity  which  is  a  turbulent  mixing 
coefficient. 

Prandtl  observed  that  v  depends  on  U  and  is  not  a 
property  of  the  fluid.  A  relation  between  v  and  the  mean 
velocity  was  needed.  Prandtl's  idea  was  that  if  small 
"lumps"  of  fluid  move  from  a  lower  to  a  higher  average 
velocity  location,  the  difference  in  its  velocity  compared 
to  the  surrounding  mean  velocity  could  be  given  by: 


where : 


U  =  f(y)  for  the  simplified  case, 

£  =  Prandtl ' s  mixing  length. 

The  physical  significance  of  l  is  the  distance  in  the  y 
direction  the  small  lump  of  fluid  must  travel  so  that  the 
difference  between  its  velocity  and  U  is  equal  to  v\-  which 
is  the  y-direction  fluctuation  velocity  at  that  location. 
After  further  development,  Prandtl  observed  that: 

-  3U  , 

-pu'v'  =  =  turbulent  shearing  stress 


This  relation  is  very  useful  in  the  calculation  of  turbulent 
flows. 

B.  GOVERNING  EQUATIONS 

The  equations  of  motion  for  the  flow  field  under 
consideration  are: 


Velocity 
field : 


V  •  v  =  0 


or 


3u  3v 
3  x  3y 


0 


Vorticity 

field: 


V  X  V  =  (A) 


or 


3  v  3u, 
3x  ‘  3y 


U) 


3u 
3 1 


+  u 


1H+Vl£=_li£  + 

3x  3y  p  3x 


v  [ 


3  2u 
3x2 


+ 


Navier-Stokes 
in  the 

x-direction : 


y^j 


which,  after  taking  the  curl  of  both  sides  of  the  equation 
and  applying  the  definitions  of  continuity  and  vorticity 
becomes  [Ref.  9]: 


=  (u-V)v  -  (vV)oi  +  vV^w 
St 


where : 


( o o*V)v  =  stretching  and  rotation 
(v-V)o)  =  convection 
vV2(jj  =  diffusion 


It  should  be  noted  that  for  turbulent  flows,  v  should  be 
replaced  by  ve,  or: 


VP  —  V  +  V 


where : 


ve  =  effective  viscosity 
vT  =  eddy  viscosity 
v  =  H  —  kinematic  viscosity. 


C.  ANALYTICAL  FORMULATION 


1 .  Kinematics 


Kinematics  is  the  branch  of  dynamics  which  deals 
with  the  motion  of  bodies  without  reference  to  the  forces 
acting  on  the  bodies.  Of  the  three  governing  equations,  the 


11 


I-,- 


n  is  the  outward  unit  normal  vector  on  B,  and 
1  1 

P  =  -  2rf  In  .+  “  which  is  the  fundamental 

I r_r0 I 

solution  for  a  two-dimensional  Poisson 
equation. 

The  consequences  of  this  formulation  include: 

The  first  integral  will  be  zero  for  the  inviscid 
region,  therefore  it  needs  to  be  computed  only  in  the 
viscous  region.  This  greatly  reduces  the 


I 


two  that  comprise  the  kinematics  of  the  flow  are  the 
equations  of  continuity  and  vorticity.  Together,  they 
express  the  relationship  between  velocity  and  vorticity 
throughout  the  fluid  at  any  point  in  time. 

There  are  two  noteworthy  aspects  of  these  relations 
[Refs.  10,11]: 

1.  The  differential  equations  are  linear,  hence  they  can 
readily  be  formulated  for  solution  by  computational 
methods . 

2.  The  stress-strain  relation  is  not  a  factor  in  these 
equations.  This  allows  the  fluid  and  solid  to  be 
treated  together  as  one  kinematic  system  which  greatly 
simplifies  the  computations. 

Through  the  use  of  fundamental  solutions  [Ref.  10] 
the  continuity  and  vorticity  field  equations  can  be 
expressed  as  [Ref.  5]: 

v(r,t)  =  -  /  w  *V  PdRQ  +  t(v0-n0)  -  (vQ  xnQ)J  *VQPdB0 

R  B 

where : 

-V 

the  subscript  "0"  refers  to  the  r0  space, 

r0  space  is  defined  by  the  vorticity  field,  or,  w0=Lo(rQ,t), 

B  is  the  boundary  of  the  region  R, 
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computational  time  required,  as  the  viscous  region  is 
only  a  fraction  of  the  total  flow  field. 

2.  The  integrals  can  be  expanded  by  Fourier  series.  This 
will  provide  more  accurate  solutions  at  each  grid 
point  than  is  possible  via  a  straight  finite 
difference  method. 

3 .  The  attached  viscous  and  detached  viscous  zone 
solutions  may  be  computed  separately.  This  allows  an 
optimized  grid  spacing  to  be  employed  in  each  zone, 
where  the  respective  length  scales  vary  greatly. 
Figure  2  presents  a  flow  zone  portrayal. 


Free  Stream  Velocity 


fnviscid  Zone  Demarcation  boundary 


Figure  2.  Flow  Zones 


4.  Repeated  computations  and  comparisons  with  either  the 
boundary  values  or  between  the  different  flow  zones 
are  not  required. 

5.  The  question  of  modeling  a  strong  viscous-inviscid 
interaction  is  obviated  completely. 

It  should  be  noted  that  all  of  these  aspects  either 
reduce  the  number  of  calculations  required,  or  enhance  the 
accuracy  of  the  solution.  The  net  result  is  an  elegant 


solution  to  the  conundrum  encountered  by  those  who  wish  to 
model  this  type  of  flow. 

2 .  Kinetics 

Kinetics  is  the  branch  of  dynamics  which  deals  with 
the  effects  of  forces  on  the  motion  of  bodies.  The 
vorticity  transport  equation  falls  into  the  realm  of 
kinetics.  It  is  nonlinear  and  elliptic  in  space.  This 
requ:  es  knowledge  of  the  solution  for  the  entire  boundary 
conditions.  For  the  outer  boundary,  this  is  met  at  the 
surface  just  inside  the  inviscid  region.  This  boundary  can 
change  with  each  time  step,  but  from  the  initial  conditions, 
the  vorticity  will  be  zero  along  this  surface.  The  interior 
boundary  is  the  solid,  and  the  vorticity  values  on  this 
surface  will  need  to  be  computed  each  time  step.  Because 
the  vorticity  values  on  the  surface  of  the  solid  are  not 
independent  of  the  vorticity  values  in  the  interior  of  the 
fluid,  this  will  have  to  be  solved  iteratively. 

3 •  Grid  Generation 

In  computational  fluid  dynamics,  grid  choice  is 
something  of  an  art.  There  are  a  number  of  conflicting 
goals  to  be  considered.  These  include: 

1.  An  adequate  number  of  grid  points  to  accurately 
represent  the  flow  conditions. 

2.  A  small  enough  number  of  grid  points  to  help  moderate 
the  computational  requirements. 

3.  Proper  grid  spacing  demands  including  a  fine  mesh  in 
the  regions  of  high  flow  variable  gradients  and  a 
coarser  mesh  in  the  other  regions.  The  latter  will 
help  moderate  the  computational  resources  required. 


4.  A  grid  that  is  easy  to  generate. 

5.  A  grid  that  simplifies  and  speeds  computations. 

Fortunately,  due  to  the  present  mathematical 
formulation,  a  very  efficient  grid  scheme  may  be  employed. 

By  choosing  the  computational  grid  to  be  circular 
with  radial  lines,  grid  generation  and  computation  goals  can 
be  met.  Because  the  flow  zones  can  be  computed  separately, 
grid  spacing  can  be  optimized  for  each  zone.  Due  to  ;he 
relatively  small  distances  from  the  solid  that  the  vortical 
flow  is  transported  to,  the  size  of  the  grid  required  can  be 
moderated.  The  only  drawback  to  this  method  is  that  the 
computational  plane  must  be  conformally  mapped  onto  the 
physical  plane.  This  limits  the  airfoil  selection  to  those 
that  can  be  accurately  mapped  between  planes.  Fortunately, 
a  number  of  airfoils  are  available  via  Joukowski  transforms, 
including  the  NACA  0012,  for  which  there  is  a  wealth  of 
previous  data  available  for  comparison. 

A  salutary  effect  of  using  a  Joukowski  transform 
between  the  computational  and  physical  planes  is  that  the 
grid  density  in  the  physical  plane  is  concentrated  radially 
about  the  leading  and  trailing  edges  while  becoming  more 
sparse  radially  over  the  upper  and  lower  surfaces  of  the 
airfoil . 

Lastly,  the  computational  grid  is  body  fixed,  which 
eliminates  having  to  recalculate  the  grid  for  each  time 
step.  Representative  grids  are  included  in  Figures  3-6. 


Initial  and  Boundary  Conditions 
Initially,  the  solid  and  the  fluid  are  at  rest.  The 
solid  is  then  impulsively  started,  and  the  discontinuity 
between  the  solid  and  the  fluid  generates  a  vorticity  sheet 
at  the  boundary.  As  time  progresses,  vorticity  is  diffused 
away  from  the  solid  and  transported  into  the  fluid  by  both 
diffusion  and  convection. 

The  four  boundary  conditions  that  must  be  met  are 
for  velocity  and  vorticity  at  the  outer  surface  of  the 
viscous  region  and  at  the  sol id/ fluid  boundary.  As  was 
previously  mentioned,  the  outer  vorticity  boundary  condition 
is  met  by  setting  the  boundary  just  inside  the  inviscid 
region.  The  outer  velocity  boundary  condition  will  simply 


be  U. 


The  inner  vorticity  boundary  condition,  as 


previously  stated,  will  need  to  be  computed  iteratively  each 
time  step,  with  the  inner  velocity  boundary  condition  found 
by  first  determining  the  tangential  velocity  component, 
which  is  due  to  the  solid's  speed,  angle  of  attack  and 
oscillatory  motion  [Ref.  5]  and  then  using  that  value  and 
the  relation  v  =  vsolid(r't)  [Ref.  13],  to  compute  the  inner 
velocity  boundary  condition. 


DESCRIPTION  OF  THE  CODE 


The  entire  runstream  has  been  named  ZETA,  which  is  an 
acronym  for  zonal  procedure  for  evaluating  turbulent  and 
laminar  flows.  ZETA  is  composed  of  three  primary  sections: 

1.  GEOM — grid  generation  and  transformation, 

2.  ZONST — main  program  in  ZETA  runstream,  and 

3.  Plotting  routines — consisting  of  pressure  computations 
and  various  plotting  options. 

The  computational  loop  in  ZONST  consists  of  [Ref.  5] : 

1.  Computing  interior  vorticity  values  by  using  the 
vorticity  and  velocity  values  from  the  previous  time 
step  in  the  vorticity  transport  equation. 

2.  Computing  new  boundary  vorticity  values. 

3.  Computing  new  velocity  values. 

Appendix  B  contains  the  version  of  the  Wu  code  used  for 
this  study  and  Appendix  C  contains  notes  on  its  employment 
at  NASA  Ames . 

A  complete  cycle  through  40  degrees  of  pitch  change  at  a 
reduced  frequency  of  .15  requires  950  time  steps  and 
approximately  800  seconds  of  CPU  time  on  the  Cray  X-MP/48 
for  an  average  of  .  842  seconds  per  time  step.  Initializa¬ 
tion  can  be  completed  with  25  time  steps.  These  are  much 
more  modest  requirements  than  codes  of  similar  ability. 


A.  GEOM 

GEOM  begins  by  reading  the  transformation  and  grid 
parameters.  The  transformation  parameters  define  the 


airfoil  that  will  be  used  while  the  grid  parameters  define 
the  size  and  stretching  coefficients  which  control  the 
radial  grid  spacing. 

GEOM  next  generates  the  computational  plane  and  then 
conformally  maps  it  onto  the  physical  plane.  The  outputs 
from  GEOM  are  used  by  all  the  other  programs.  Grid  point 
spacing  may  be  checked  here  and  adjusted  as  necessary. 


B.  ZONST 


ZONST 

is  the  main 

program  of 

ZETA. 

It 

uses  the 

governing 

equations  to 

compute  the 

vorticity 

and 

velocity 

fields  and  generates 

the  output 

used  by 

the 

plotting 

routines . 

It  starts  by  reading  the  grid  information  from  GEOM  and 
its  own  input  parameters.  The  input  options  available 
include: 

1.  Airfoil  motion.  This  can  be  defined  as  a  constant 
angle  of  attack,  a  rapidly  pitching  or  an  oscillating 
motion.  All  pitching  airfoil  flow  solutions  require 
data  from  the  appropriate  constant  angle  of  attack 
steady  state  case  to  be  stored  as  part  of  the  initial 
conditions . 

2.  Time  specifications  allow  the  user  to  control  the 
number  of  time  steps  to  be  run;  the  increment  of  the 
time  steps  and  the  time  step  values  fur  which 
numerical  and  graphical  output  will  be  generated. 

3.  Flow  zone  specifications  can  be  used  to  control  the 
boundary  layer  region  where  Navier-Stokes  computations 
will  be  used.  Generally,  there  will  be  no  increase  in 
accuracy  but  a  significant  increase  in  computation 
time  when  using  Navier-Stokes  vice  boundary  layer 
calculations  in  the  boundary  layer  region. 


20 


4.  Under-relaxation  parameters  allow  tailoring  or  the 
computational  methods  to  help  optimize  computing 
efficiency  for  various  flow  conditions. 

20NST  can  be  run  for  a  specified  period,  checked  and 
restarted.  Should  divergence  occur,  the  last  computed  time 
step  will  be  printed  for  trouble  shooting. 

The  output  from  ZONST  is  numerical,  and  via  plotting 
routines,  graphical. 

C .  PLOTTING  ROUTINES 

Due  to  the  primary  variables  used  in  ZONST,  the 
aerodynamic  loads  caused  by  the  pressure  distribution  are 
not  directly  available.  LOADS  uses  the  vorticity 

information  and  Fourier  series  expansions  of  the  vorticity 
integral  representation  to  derive  values  for  coefficients  of 
pressure,  lift,  drag  and  moment  for  each  angle  of  attack. 

The  first  three  plotting  routines  listed  use  DISSPLA  for 
their  operating  software,  while  the  last  one  uses  PL0T3D . 

Plotl  portrays  the  physical  plane  vorticity  grid,  the 
boundary  layer  grid,  flow  zone  demarcations  and  wake  and 
turbulence  grids. 

Plot2  generates  streamline  and  vorticity  contours  as 
well  as  numerical  output  of  the  contours  and  grid  points 
crossed . 

Plot3  is  the  loads  plotting  program.  Coefficient  of 
pressure  versus  non-dimensionalized  chord  length  for  each 
selected  angle  of  attack  may  be  displayed,  or  plots  for  the 
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RESULTS  AND  DISCUSSION 


A.  COMPARISON  WITH  EXPERIMENTAL  DATA 

The  first  comparison  was  with  experimental  data  from 
[Ref.  14]  for  a  series  of  three  different  reduced 
frequencies,  at  the  same  Mach  and  Reynolds  numbers.  This 
was  selected  to  highlight  the  time  history  dependent  nature 
of  dynamic  stall.  The  plots  pertaining  to  this  comparison 
are  enclosed  in  Appendix  A.  Considered  were  Cp  versus  non- 
dimensionalized  chordlength  or  x/c,  as  well  as  CL,  Cm  and  Cp 
versus  incidence  angle,  a. 

Reduced  frequency  is  a  pitching  rate  parameter  and  is 
defined  as: 


_  COC 

rf  =  2C 


where : 


co  =  circular  frequency, 

C  =  airfoil  chord  length,  and 
U  =  freestream  velocity. 


The  experimental  data  were  taken  at  M  =  .072  to  provide 
a  valid  approximation  of  incompressible  flow. 


For  rf  =  .099,  the  plots  of  Cl,  Cp  and  Cm  versus  ct  pro¬ 


vide  a  good  overview  of  the  conditions  the  airfoil 


experiences.  In  all  three  plots,  the  theoretical  data 
follows  the  trends  of  the  experimental  data  well,  while  the 
relative  magnitudes  and  short  term  stall  recovery  features 
are  less  well  represented. 

The  theoretical  values  of  Cp  also  track  relatively  well 
with  the  experimental  data,  though  slightly  under-represent¬ 
ing  the  pressures  recorded  experimentally.  The  theoretical 
results  indicate  the  airfoil  as  beginning  to  stall  at  a 
slightly  lower  angle  of  attack  than  the  experimental  results 
do.  At  01  =  20.8,  the  theoretical  data  show  the  vortex 

bubble  as  having  been  shed  from  the  leading  edge  of  the 
airfoil.  Its  progression  across  the  upper  surface  can  be 
followed  in  subsequent  plots. 

Streamline  and  vorticity  contour  plots  are  presented  to 
illustrate  the  physical  phenomena  but  are  not  directly 
compared  with  experimental  data.  Vorticity  bubble 

generation  and  propagation  are  observed  in  each  of  the  three 
cases . 

For  the  plots  of  rf  =  .15,  stall  onset  is  noticably 

delayed  from  the  rf  =  .099  case.  Stall  now  occurs  at  a  ' 
24°  for  both  the  experimental  and  theoretical  cases.  As  in 
the  rf  =  .099  case,  the  trends  show  good  correlation  with 

post  stall  effects  being  less  well  modeled. 

For  rf  =  .248,  the  plots  of  C l,  Cq  and  Cm  versus  a  have 
somewhat  poorer  correlation  than  the  previous  two  cases. 


There  is  greater  offset  between  data,  but  stall  occurs  at 
a  =  25°  for  both  data  sets  on  the  CL  versus  curves. 

For  these  three  cases,  there  is  a  slight  trend  towards 
higher  maximum  values  on  the  CL  versus  ot  curves  with  higher 
reduced  frequency.  Significantly,  the  theoretical  results 
accurately  reflect  the  trend  towards  a  higher  stall 
incidence  angle  with  higher  reduced  frequency. 


B.  COMPARISON  WITH  OTHER  THEORETICAL  DATA 

For  the  last  case,  a  comparison  was  made  between 
experimental  data  and  two  types  of  theoretical  data. 
Experimental  data  are  from  [Ref.  14]  and  one  version  of  the 
theoretical  data  is  from  [Ref.  15].  The  conditions  for  this 
case  were:  reduced  frequency  =  .199,  mach  number  =  .184  and 
-2  <  a  <  18°. 

Figure  7  is  the  experimental  and  other  theoretical  data. 
Figure  8  is  the  theoretical  data  from  the  Wu  code. 
Referencing  the  Cp  versus  x/c  plots  in  Figs.  7a, b  and  8c, 
the  theoretical  data  from  the  Wu  code  demonstrate  a  very 
high  degree  of  accuracy  with  the  experimental  data,  while 
the  theoretical  data  from  the  other  source  show  a  decidedly 
less  accurate  picture.  The  other  source  indicates  dynamic 
stall  and  post  stall  conditions  when  experimentally,  the 
airfoil  never  stalled. 

Similar  and  pronounced  differences  also  occur  in  the  CL 
and  Cm  versus  a  plots.  As  in  the  previous  series,  the  Wu 
code  data  slightly  under-represents  the  amplitude  of  the 
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pressure  but  still  follows  the  trends  well,  while  the  other 
theoretical  data  which  incorrectly  indicates  a  dynamic  stall 
condition. 

For  this  case,  the  data  from  the  Wu  code  data  are  a  more 
accurate  reflection  of  the  experimental  data  than  the  data 
from  the  other  theoretical  source. 


C.  ENHANCEMENTS  TO  THE  CODE 

There  were  two  primary  enhancements  that  were  added  to 
the  Wu  code.  These  are: 

1.  The  ability  to  output  plotting  data  at  non-regular 
intervals,  to  add  flexibility  to  the  code  and  to 
simplify  and  facilitate  comparison  with  experimental 
data. 


2.  The  availability  of  the  code  to  generate  physical 
plane  velocity  data  to  allow  the  use  of  Plot3D  and 
its  associated  graphics  functions. 

The  modifications  to  provide  physical  plane  velocity 
information  centered  around  the  transformation  matrix  that 
is  used  in  GEOM  to  obtain  the  physical  plane  grid  from  the 
computational  plane  grid  after  the  computational  grid  is 
generated. 

The  scale  factor  of  the  transformation  is  defined  as 


where  Z  is  the  physical  plane  and 


V 
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while  z,  is  the  computational  plane,  where  the  airfoil  is 
represented  by  the  unit  circle  and  cylindrical  coordinates 
are  employed. 

$  z 

Then,  denoting  as  DZDW,  HSTAR  is  defined  as: 

H*  =  complex (real (DZDW) , -absolute  imaginary  (DZDW)) 

HSTAR  is  written  onto  tape  20  and  passed  to  ZONST,  where 
it  is  used  in  the  subroutine  KMTCS.  In  KMTCS,  two 
additional  arrays  are  defined,  VTOTCO  and  VTOTPH,  which  are 
the  total  computational  plane  and  physical  plane  velocities, 
respectively.  They  are  defined  as: 

VTOTCO  =  complex  (W1,W2)  and 
VTOTPH  =  VTOTCO/HSTAR 

where  Wl  and  W2  are  the  computational  plane  velocities  in 
the  rotating  frame  of  reference.  Additional  minor  modifica¬ 
tions  to  the  code  were  made  to  implement  these  changes. 

Representative  plots  are  presented  for  four  cases  where 
the  reduced  frequency  is  .150  and  Reynolds  number  is  1  x106. 

1.  Laminar  flow  with  the  airfoil  at  a  steady  state  5.0° 
a.  This  is  portrayed  in  Fig.  9. 


IS: 


2.  Flow  conditions  immediately  prior  to  stall  onset  are 
shown  in  Fig.  10. 

3.  Initial  indication  of  reverse  flow  at  21.06°  i  in  Fig. 
11.  This  appears  on  the  upper  surface  near  the 
leading  edge  as  the  vorticity  bubble  starts  to  form, 
which  quickly  spreads  over  the  entire  upper  surface  of 
the  airfoil. 

4.  Turbulent  flow  with  the  airfoil  at  23.8°  a,  shortly 
after  the  onset  of  dynamic  stall.  This  is  portrayed 
in  Fig.  12.  In  all  cases,  alternate  radial  grid  lines 
were  deleted  from  the  plots  to  enhance  clarity,  while 
a  similar  display  progression  is  followed  for  the 
first  and  last  cases. 

Two  additional  modifications  would  help  interpretation 
of  the  information.  First,  computing  VTOTPH  for  the  entire 
physical  plane,  not  just  in  the  viscous  region.  Second, 
rotating  the  grid  so  that  the  airfoil  is  graphically  shown 
at  its  true  angle  of  incidence.  Both  of  these  deficiencies 
are  presented  in  Fig.  9a,  and  neither  affects  the  accuracy 
of  the  graphical  output  in  the  viscous  region. 

In  Fig.  9a,  the  primary  portion  of  the  viscous  region 
and  part  of  the  inviscid  region  are  shown. 

Subsequent  portions  of  Fig.  9  show  increasing  resolution 
of  various  portions  of  the  plot  in  Fig.  9a.  To  aid  in 
comparison,  the  horizontal  scale  is  consistent  throughout 
the  velocity  plots. 

Fig.  9c  shows  the  leading  edge  stagnation  point  and 
tangency  of  the  boundary  layer  at  the  surface  of  the 
airfoil.  Fig.  9d  is  the  midchord  upper  surface  region. 
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detail  plot  of  the  boundary  layer  near  the  mid-span,  upper 


surface . 


Fig.  11  is  the  case  for  dynamic  stall.  Flow  reversal  is 
apparent  over  the  entire  upper  surface,  with  the  vortex 
bubble  center  indicated  by  the  zero  velocity  vector  located 
above  approximately  .3  chord.  In  Fig.  11c,  the  high 
velocity  gradients  are  readily  apparent,  from  the  leading 
edge  around  to  the  upper  surface. 

More  complex  flow  patterns  can  be  readily  portrayed. 


V.  CONCLUDING  REMARKS 


The  Wu  code  holds  much  promise  to  help  unlock  some  of 
the  mysteries  of  dynamic  stall.  Its  strengths  include  the 
following: 

-  Enough  speed  and  efficiency  that  it  can  be  operated  on  a 
VAX  or  similarly-sized  computer. 

-  Good  success  at  indicating  the  trends  of  the  cases 
studied. 

-  Relatively  accurate  results. 

-  Powerful  diagnostic  tool  which  can  become  a  predictive 
tool . 

-  Not  Reynolds  number  limited. 

Other  aspects  that  should  be  noted  are: 

-  The  mathematical  formulation  is  more  involved  than  a 
straight  finite  differencing  of  pressure/velocity 
Navier-Stokes  equations. 

-  As  currently  formulated,  it  is  not  readily  applicable  to 
arbitrary  geometries,  but  a  useful  selection  of 
Joukowski  transforms  is  available. 

With  the  addition  of  compressibility  and  transition 
modeling  and  later,  extension  to  three  dimensional 
representation,  its  utility  will  continue  to  be  enhanced  and 
its  applications  expand. 


APPENDIX  A 


REDUCED  FREQUENCY  VARIATION  PLOTS 


Appendix  A  contains  plots  for  variation  of  reduced 
frequency.  For  each  of  the  three  reduced  frequencies,  the 
order  is:  C £,  CD  and  Cm  versus  a 

Cp  versus  x/c  with  varying  a 

Streamline  and  vorticity  contours  with  varying  a. 
Theoretical : 


Experimental : 
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APPENDIX  B 


CODE  LISTING 


PROGRAM  GEOM 


GEOMETRY  PLOTTING  PROGRAM  -  DISSPLA  VERSION 
LAST  REVISION  4-1-87 

PRINCIPAL  INVESTIGATOR  OD  DR.  J.C.  MJ 
AUTHORS  00  MIKE  PATTERSON.  ISHMAEL  TUNCER 
GEORGIA  INSTITUTE  OP  TECHNOLOGY 
(404)  894-3028 


TAPE1 

OD 

INPUT  TO  PL0T1 

TAPE2 

00 

input  TO  ZONST 

TAPE3 

OD 

INPUT  TO  LOADS 

TAPE5 

00 

GENERAL  INPUT 

TAPE6 

OD 

GENERAL  OUTPUT 

TAPE10 

OD 

XYZ  INPUT  TO  PL0T3D 

TAPE14 

00 

INPUT  TO  PL0T2.  PLOT 3 

CALLS 

00 

«#*< 

NONE 

IMPLICIT  REAL'S  (A-H.O-Z) 

PARAMETER  ( I DIM-80 , JDIM-60) 

PARAMETER  ( I PI-81 . JP 1-61 ) 

PARAMETER  (KFC1-41 .KFC2-42) 

PARAMETER  (INOR1-20) 

DIMENSION  UC(IP1 ,JP1) .VC(IP1 . JP1 ) .HREAL(IDIM) .HIMAG(JDIM) 
DIMENSION  CS(KFC1 . ID1M) ,SN(KFC1 . IDIM) ,HSTAR( IDIM, JDIM) 
DIMENSION  R1 (JP1 ) ,R2( JP1 ) ,R1D( JP1 ) ,RP( JP1 .KFC2) ,RL(JP1 ) 
DIMENSION  H(IDIM.JDIM) .SGMA(KFCI) 

DIMENSION  A2(JDIM) ,A4( JDIM) ,C2( JDIM) .C4( JDIM) , D5( JDIM) 
DIMENSION  AS2(KFC1).BS2(KFC1) ,CS2(KFC1) ,DS2(KFC1) 
DIMENSION  YN(KFC1 , JDIM) .CCP(KFC1 .JDIM) 

DIMENSION  I WK( JDIM, JDIM) . INOR( IN0R1 , JDIM) ,COEF( IDIM. 2) 
DIMENSION  X( IDIM. JP1 ) ,Y( IDIM, JP1) .XB( IDIM) , YB( IDIM) 
DIMENSION  XP(IP1 . JP1 ) . YP(IP1  , JP1 ) ,XT( IDIM) , YT( IDIM) 
COMPLEX  W.*1 .Z.Z1 .DZDW.OWT.HIAMG.HSTAR 

READ(5 , •)  RE 

READ! 5 . •)  C1.CC2.DS 

READ ( 5 , • )  CSO.GAMA.SIGMA.AL 

IM-IDIM 

JR— JDIM 

JR1-JP1 

KFC-KFC1 

PI— 4 . »ATAN( 1 .) 

OTET— 2. »PI/FLOAT(IM) 

I M2- I M/2 

C.. COMPUTE  RAOIAL  GRID  DISTRIBUTION  IN  COMPUTATIONAL  PLANE 

00  100  J-1 . JR 1 
S-FLOAT( J-1 )»0S 


R1 ( J )«EXP(S+C1 )+CC2 
R2(J)»EXP(S+C1+.5«DS)+CC2 
100  CONTINUE 

DO  105  J-1 .JR 
105  R1D(J)-R1 (J+1 )-R1 ( J) 

C.. COMPUTE  SCALE  FACTOR  OF  TRANSFORMATION 

DO  110  J-1 , JR 
DO  110  1-1 . IM 

PHI«FLOAT( 1—1 )*DT£T 
WR— R2( J )»COS(PHI ) 

WI-R2(J)*SIN(PHI) 

W-CMPLX(WR.WI) 

ZI-WlCAMA 

DZDHM  .-CSQ/(Z1  *Z1 ) 

HI-OABS(DZDW) 

H(  I ,  J)— HI  »H1 
HREAL(I)  -  REAL(DZDW) 

HIMAG(J)  -  AIMAG(DZDW) 

HSTAR(I.J)  -  CMPLX(HREAL( I ) ,-HIMAG( J) ) 

110  CONTINUE 

WRITE(20)  HSTAR 

C.. ARRAYS  USED  IN  KINETICS 
DO  115  J-1 .JR 

A2( J)-R1 (J+1 )/(DS*(R2( J )-CC2) ) 

A4( J )— A2(  J ) *AL/(RE»DS* (R1 (J+1 )-CC2) ) 

C2(J)— R1(J)/(DS*(R2(J)-CC2)) 
C4(J)-C2(J)*AL/(RE*DS*(R1(J)-CC2)) 

D5( J )-AL/(R2( J ) •RE«DTET »OTET } 

115  CONTINUE 

C. .COMPUTE  GEOMETRIC  COEFFICIENTS  FOR  COMPUTATION  OF  VELOCITY 

C  FOURIER  COEFFICIENTS  BY  INTEGRAL  RELATIONS 

DO  120  J-1 , JR1 
RP(J , 1 )— R1 ( J) 

RL(J)-LOG(R1(J)) 

DO  120  K— 2.KFC+1 
RP(J.K)— RP(J,K-1).RP(J.1) 

120  CONTINUE 

C.. COMPUTE  COORDINATES  AND  DERIVATIVES  OF  SOLID  SURFACE 

C  FOR  USE  IN  LOADS  CALCULATIONS 

DO  125  1-1 .IM 

PHI«FLOAT( I— 1 )»DTET 
CO-COS(PHI) 

SI-SIN(PHI) 

W-CMPLX(CO.SI) 

DWDT-CMPLX(-SI  .CO) 

Z1-W+GAMA 

Z-Z 1 +CSQ/Z 1 +S I GMA 

XB(I)-REAL(Z) 

YB( I )— AIMAG(Z) 

DZD»N1 ,-CSQ/(Z1 *Z1 ) 

XT(I)-REAL(DZDW.DWOT) 

YT(  I )-AIMAG(DZDW*DWDT) 

125  CONTINUE 

C.. COMPUTE  CORRELATIONS  OF  VELOCITY  COMPONENTS  BETWEEN 

C  INERTIA  COORDINATE  SYSTEM  AND  BODY  COORDINATE  SYSTEM 


00  133  1=1 . IM 

PHI-FLOAT ( 1-1 )»DTET 
WR«R1(J).C0S(PHI) 

V*I-«1(J).SIN(PHI) 

W1-WR+GAMA 

W2-W1-W14WI-W 

XR-W1+CSQ*W1/W2+SIGMA 

YR-WI-CS0*WI/X2 

XOI  .+CSQ»(WI*WI-W1«W1)/(W2*W2) 

XN—CSQ*  2.«Wl«W1/  (W2*W2  ) 

UC(I .J)-XR»XN+YR*XC 
VC(I.J)-YR*XN-XR*XC 
130  CONTINUE 

DO  135  J-1 , JR1 

UC( IM41 , J)— UC( 1 , J) 

VC(IM4-1  ,J)-VC(1  ,J) 

135  CONTINUE 

C.. COMPUTE  COORDINATES  OF  VORTICITY  GRID  IN  PHYSICAL  PLANE 
C  AND  IN  COMPUTATIONAL  PLANE 


DO  140  J-1 , JR1 
DO  140  1—1 . IM 

PH I -FLOAT (1—1 )*DTET 
WR— R2( J) *COS(PHI ) 

WI-R2(J)*SIN(PHI) 

W-CMPLX(WR.WI) 

Z-#+GAMA4CSQ/(W+GAMA)+SIGMA 

XP( I , J)— WR 

YP(I.j)-WI 

X(I.J)-REAL(Z) 

y(i  ,  j)— aimag(z) 

140  CONTINUE 

DO  141  J-1.JR1 
XPfIMM-l . J)— XP(1 ,  J) 

YP(  IU+-1  .  J  )  -YP  ( 1  .J) 

141  CONTINUE 
WRITE(10)  IDIM.JDIM 

WRITE( 10)  ( (X( I . J) . 1-1 . IDIM) . J-1 .  JDIM) , 
1  ((Y(I.J).  1-1. IDIM), J-1, JDIM) 

DO  145  1-1 .IM 

PHI— FLOAT(I— 1 )*DTET 
CS(1.I)-1. 

SN(  1  . 1 )— 0 
CS(2.I)-C0S(PHI) 

SN(2 . I )— SIN(PHI ) 

145  CONTINUE 

DO  150  K-3.KFC 
DO  150  1-1 .IM 

L-14M0D( (K-1 ) • ( 1-1 ) . IM) 
CS(K.I)-CS(2.L) 

3N(K . I )— SN(2 . L) 

150  CONTINUE 

C.. COMPUTE  FOURIER  SERIES  SMOOTHING  FUNCTION 


DO  155  K-2.KFC 

A1«FL0AT(K-1)»0TET 
SCMA(K)— SIN(A1 )/A1 
155  CONTINUE 


SWij 
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C.. DETERMINE  FOURIER  COEFFS.  OF  BOUNDARY  VORTICITY  CONDITION 
C  ON  INERTIA  COOR .  (COEF(I.I)  FOR  RADIAL  COMP.,  COEF(l,2)  FOR 
C.  TANGENTIAL  COMP.) 


DO  160  K-1 , KFC 
AS2(K)=0. 
BS2(K)-0. 
CS2(K)-0. 
DS2(K)-0. 

160  CONTINUE 

DO  161  1-1. IM 


COEF( I , 1 )— (UC( I . 1 )*CS(2, I )+VC( I . 1 )*SN(2 , I ) )  /  FLOAT ( I M2) 
COEF( I , 2)— (VC( I , 1 )»CS(2, I )-UC( I , 1 )»SN(2, I ) )  /  FLOAT ( I M2) 


161  CONTINUE 

DO  162  1-1. IM 

AS2( 1 )— AS2( 1 )+COEF( 1.1) 
CS2(  1 )-CS2( 1 )+CO£F( 1,2) 


DO  162  K-2.KFC 

AS2(K)— AS2(K)+COEF ( 1.1) «CS(K , I ) 
BS2(K)-8S2(K)+C0EF( I . 1 )»SN(K , I ) 
CS2(K)«CS2(K)+CO£F ( I , 2) *CS(K . I ) 
DS2(K)-OS2(K)+CO£F ( I , 2) *SN(K , I ) 

162  CONTINUE 

DO  163  K-2.KFC 

AS2(K)— AS2(K) *SGMA(K) 
BS2(K)=8S2(K)*SGMA(K) 
CS2(K)-CS2(K)«SGMA(K) 

DS2(K)« 0S2(K)*SGMA(K) 

163  CONTINUE 


C. .COMPUTE  ‘ YN‘ 


DO  165  1-1. KFC 
DO  165  J-1.JR 
IS-I 

I F( I  .GT.  I M/3)  THEN 

YN( I . J)— SQRT( (X( I , J )-XB( I ) ) ••2+{Y( I , J )— YB( I ) ) *»2) 
ELSE 

0 1 F— X  ( I  .  J  )— XB  ( I  ) 

I F ( ABS (DIF)  .GT.  0.001)  THEN 
INC-1 

IF(DIF  .GT.  0.)  INC— 1 

164  IS-IS+INC 

I F( IS  GT.  KFC  .OR.  IS  -LT.  1)  THEN 
YN(I.J)-Y(I,J) 

ELSE 

DIFN-X(I,J)-XB(IS) 

IF(DIFN*OIF  -LT.  0.)  THEN 
ISU-IS-INC 

YIS— YB( ISM)+(X( I , J )-XB( ISM) )• (YB( IS)— YB( ISM) ) 
1  /(XB(IS)-XB(ISM)) 

YN(I.J)— Y(I ,J)-YIS 
ELSE 

GO  TO  164 
ENDIF 
END  IF 
ELSE 

YN( I , J )»Y( I , J )-YB( I ) 

ENDIF 

ENDIF 

165  CONTINUE 


C. .COMPUTE  1 INOR’ 


"J*  V  V*  V  ", 


,"v  ivmvv 


DO  170  II-2.IM/4 
I-I I 

INOR( 11,1 )— I 
DO  170  J-2.JR 
INOI 

I F(X( I ,  J )  .LT.  X(1 1 . 1 ))  INC— 1 

IF(A8S(X( I 1 , 1 )— X( I+INC, J) )  .LT.  ABS(X( 1 1 . 1 )-X( I . J ) ) )  I-I+INC 
INOR(II ,J)«I 
170  CONTINUE 

C.. COMPUTE  ’INK’ 

DO  175  JKt-1  ,  JR— 1 
I-IN0R(2.JN) 

IWK(JN.JN)-I 
DO  175  J-JN+1.JR 
INOI 

IF(Y(I,J)  .GT.  Y(I.JN))  INC— 1 

IF(ABS(Y(I.JN)-Y(I+INC,J))  .LT.  ABS(Y( I . JN)-Y( I . J) ) )  I-I+INC 
IWK(JN, J)-I 
175  CONTINUE 

C.. CALCULATE  GEOMETRIC  COEfT.  ON  PRESSURE  COMP. 

DO  180  J-1 .JR 
DO  179  I-3.KFC 

179  CCP(I.J)-  I-  /  FLOAT(I-2).(1./R1(J)**(I-2) 

1  — 1 ./R1 ( J+1 )**(I-2) ) 

CCP(1 .J)-  R1D(J) 

CCP(2.J)-  LOG(R1(J+1)/R1(J)) 

180  CONTINUE 

C.  INTEGRATE  OVER  AIRFOIL  SURFACE  TO  COMPUTE  AIRFOIL  AREA; 

C  ASSUMES  AIRFOIL  IS  SYMMETRIC 

DPHI-PI/200. 

CO— 1. 

SI-0. 

*-<MPLX(CO.SI) 

Z1-W+GAMA 

Z-Z1+CSQ/Z1 

XI-REAL(Z) 

Y1— AIMAG(Z) 

PHI-PI 
A1— 0. 

DO  200  1-1 .200 
PHI— PHI— DPHI 
CO-COS (PHI) 

SI-SIN(PHI) 

W-CMPLX(CO.SI) 

Z1-N+GAMA 
Z-Z1+CSQ/Z1 
X2— REAL(Z) 

Y2— AIMAG(Z) 

A1-A1+.5*(Y2+Y1)«(X2-X1) 

XI -X2 
Y1-Y2 

200  CONTINUE 
AF— 2 • *A1 


C.  OUTPUT  TO  THE  VARIOUS  FILES 
NRITE(I)  AL.X.Y, INOR . INK 


w 


k 


WRITE(2)  DTET.DS.AF.AL.RE 
WRITE(2)  CS , SN , UC , VC , H 
WRITE (2)  RP.RL.SGMA 
WRITE(2)  A2.A4,C2,C4,D5 
WRITE(2)  AS2 , BS2 , CS2 , DS2 
WRIT£(2)  R1 ,R2,R1D, YN 
WRITE(2)  I NOR , IWK ,CCP 

WRITE(3)  AL , RE , DTET , XB , YB , XT , YT 
WRITE(3)  CSQ.GAMA, SIGMA, Cl ,CC2 

WRITE(6.10) 

WRITE(6,40)  IM, JR, SIGMA, CSQ. Cl .CC2.DS.GAMA.AL.RE 
WRITE(6,20)  (H( 1.1) ,1-1 ,IM) 

WRITE(6.30)  R1 

WRITEM4)  CSQ.GAMA.SIGMA, AL 
WRITE(14)  XB.YB.XP.YP 

10  FORMAT (//, ’INPUT  TO  GEOM’) 

20  FORMAT (//, 'SQUARE  OF  THE  SCALE  FACTORS  AT  FIRST  RINGOO’ 

1  ,//( 10E13 . 6) ) 

30  FORMAT (// , ’  GRID  DISTRIBUTION  IN  R  DIRECTION  OD  ’//(10E13.6)) 

40  FORMAT (//./IX, 

1  ' IM«  ’ .T10.I5./1X, 

2  'JR-  ’ .T10.I5./1X, 

3  'SIGMA-  ' .T10.F10.7./1X, 

4  'CSQ-  '  .T10.F10. 7,/IX. 

5  'Cl-  ' , T10 , F10. 7,/IX , 

6  'CC2—  ' .T10.F10.7./1X. 

7  'DS-  ' .T10.F10.7./1X, 

8  'GAMA-  ' .T10.F10. 7,/IX, 

9  'AL-  ' .T10.F10.7./1X, 

1  'RE-  '  . T10 , F10. 1 ) 

STOP 

END 
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1 000000 

-4.  9279788*7591  9927588760838  .1293603008237 

.8061367  -.0647384  .9170954997604  3.619058974007 


INPUTS: 


NACA  0012  REVISED  PARAMETERS 


RE 

Cl .C2.DS 

CSQ.GAMA.S1GMA.AL 


u 


000000000000000000000)0 


PROGRAM  ZONST 


2-0  INCOMPRESSIBLE  NAVIER-STOKES  SOLVER  -  DISSPLA  VERSION 
LAST  REVI SIONOO  5-20-87 

PRINCIPAL  INVESTIGATOR  00  OR.  J.C.  WU 
AUTHORS  00  MIKE  PATTERSON,  ISHMAEL  TUNCER 
GEORGIA  INSTITUTE  OP  TECHNOLOGY 
(404)  894-3028 


TAPE2 

00 

OUTPUT 

FROM  GEOM 

TAPE4 

00 

OUTPUT 

FOR  LOADS 

TAPES 

00 

GENERAL 

INPUT 

TAPE6 

00 

GENERAL 

OUTPUT 

TAPE7 

00 

INPUT  FROM  PREVIOUS 

TAPE8 

00 

OUTPUT 

FOR  NEXT  RUN 

TAPE9 

OD 

OUTPUT 

FOR  PLOT 2 

TAPE11 

00 

Q  INPUT 

FOR  PL0T3D 

CALLS 

00 

KMTCS 

KNTCS 

CPVAL 

PARAMETER  ( IDUA-80,  JDIM-60) 

PARAMETER  ( IP1-81 . JP1-61  ) 

PARAMETER  (KFC1-41 .KFC2-42) 

PARAMETER  ( IDUfc^4295 . INOR1-20) 

DIMENSION  AS2(KFC1 ) .BS2(KFC1 ) .CS2(KFC1 ) ,DS2(KFC1 ) 
C0M40N/I0/L00P , NT , NTMAX , NTOUT 

C0M40N/DKIN/A2( JDIM) ,A4(JDIM) ,C2( JDIM) ,C4(J0IM) ,D5( JDIM) 
COfcMON/V  LB/VOR  LB ,  NLB 
COMMON/RGRD/R1D( JP1 ) ,R2(JP1 ) 

C0MM0N/SMT/SGMA(KFC1 ) 

COMMON/RFC/AA( JDIM, KFC1 ), BB( JDIM, KFC1 ) 

C0M40N/RHS/#0W(  IP1  .  JP1  )  ,POP(  IP1  .  JPl  )  ,WP(  IP1  )  ,WB( IP1  )  ,R1  (  JP1  )  , 

1  DUMM(IOUM) 

C0M40N/ ABCDS/AS 1 (KFC1 ) ,BS1 (KFC1 ) fCS1 (KFC1 ) ,DS1 (KFC1 ) 
COMMON/COR/RP ( JP 1 ,KFC2) ,RL( JPl ) 

C0M40N/0ELTA/DS , DTET , OT 
COfcMON/UNF/URR . URF . DFMX , NCC 
C0»*40N/VTEXT/KIN(I0IM)  .KEN(IOIM) 

C0MM0N/TRIG/CS(KFC1 . IDIM) , SN(KFC1 , IDIM) 

C0M40N/SCALE/H( IDIM , JDIM) 

C0MM0N/W/V0R(  IDIM,  JDIM)  ,VOROLD(  IDIM.  JDIM) 

C0M40N/VEL/U(IDIM, JPl), V( IDIM, JPl ) , HSTAR( IDIM. JDIM) 
C0MM0N/TUR1/USTAR( IDIM), EDOY( IDIM, JDIM) . IOTUR(IDIM) . 

1  YN(KFC1 . JDIM) , IOT. IOTB 

C0MM0N/TUR2/ 1 NOR ( I  NOR 1 .JDIM) . IWK( JDIM. JDIM) 

COMMON/VEH/UC ( IP1 .JPl ) ,VC( IP1 , JPl ) 

C0M40N/ZNS/IB1 .IV1 . IV2.IB2.IV1R 
COMMON/GRD/IM, IM2.KFC.JR.N 

C0M40N/C0E/AF.UI .VI .OMG . VSC . NPL . ICTUR . ICST , ICPL 
COMMON/CPC/CCP (KFC 1 . JD IM) , CP ( I P 1 ) 

C0M40N  Q1 (IDIM, JDIM) ,Q2( IDIM, JDIM) ,Q3( IDIM. JDIM) ,Q4( IDIM. JDIM) 
COMPLEX  HSTAR 

REAL  ALLP(2) 

DATA  ALLP/  19.5.20.0/ 
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V 


tv* 
$ 
$ 


.-I 


f.  -t 


km 


'.v, 


ifijl 


Si? 


:Si 

A 


o  o  o  o  o  o  o 


REWIND  2 
REWIND  4 
REWIND  5 
REWIND  6 
REWIND  7 
REWIND  8 
REWIND  9 
REWIND  20 

C. .READ  GENERAL  INPUTS  AND  ECHO  THEM 

READ(5, •}  ICST.RF.ALPS. ICTUR 
READ(5, • )  WMIN,DFMX,DRMX,KMAX,NCC 
READ(5, •)  URBI ,URBP , URR 
READ(5, • )  IV1 , IB1 , IB2, IV2,NPL,NLB 
READ(5, •)  DTI .DTINC.NTMAX 
READ(5, •)  NTPL.NTOUT.NTLO 

WRITE(6.54) 

WRITE(6, •)  ICST.RF.ALPS, ICTUR 
WRITE(6 . • )  WMIN.DFMX ,DRMX , KMAX , NCC 
WRITE(6, *)  URBI, URBP. URR 
WRITE(6. •)  IV1 ,IB1 . IB2, IV2.NPL.NLB 
WRITE(6 . • )  DTI .DTINC.NTMAX 
WRITE(6, •)  NTPL.NTOUT.NTLO 


C. . INPUTS  FROM  GEOM 

READ(2)  DTET.DS.AF.AL.RE 
P?EAD(2)  CS.SN.UC.VC.H 
READ(2)  RP.RL.SGMA 
REAO(2)  A2.A4.C2.C4.05 
READ(2)  AS2.BS2.CS2.DS2 
READ(2)  R1.R2.R10.YN 
READ(2)  INOR.IWK.CCP 
READ(20)  HSTAR 

NALLP-1 

IM-IDIM 

JR-JDIM 

KFC-KFC1 

URF-1  -URR 

NLBP-NLB+1 

PI-4. *ATAN( 1 . ) 

IM1-IU+1 

IV1R-IM4-IV1 

IM2-IM/2 

RDTET-DTET.R1 (NLBP) 

VSC-AL/RE 
FR-2. *RF/AL 
ALP-ALPS »P I /1 80. 

OMG-0. 

OMGD-0. 

C.. START  INITIAL  SOLUTION  OR  READ  PREVIOUS  ITERATION  RESULTS 

IF( ICST . LT . 2)  THEN 
DO  100  K-1 . KFC 
AS1  (K)— 0 . 

BS1 (K)-O. 

CS1 (K)» 0 . 

DS1(K)-0. 

100  CONTINUE 


END  IF 

I F ( ICST . EQ. 0)  THEN 
NT-0 

I F(DT I .EQ.0. )  DT-0.08 

T-0.0 

N— 1 5 

VORLB-0 . 0 
UI-COS(ALP) 

VI-SIN(ALP) 

00  105  1-1 . IM 
KIN(I)«N 
U(I.1)«0. 

V(I .1)-0. 

DO  105  J-1 , JR 
VOR( I , J )— 0 . 

105  CONTINUE 
AA(1 , 1 1-0. 

AA( 1 ,2)— 2. *Vl/(RP(2, 1 )— RP(1 , 1 )) 

BB(1 ,2)— 2. *UI/(RP(2, 1 )-RP( 1.1)) 

C.. POTENTIAL  FLOW  SOLUTION 

CALL  KMTCS 
DO  106  1-1 . IM 

106  VOR( I . 1 )«(AA(1 ,2) »CS(2 , I )+BB( 1  .2) *SN(2 . I ) )/H( I ,  1  ) 

END  IF 

C  .READ  PREVIOUS  RUN  AND  RESET  TIME  IF  A  NEW  MOTION  IS  SPECIFIED 
I-ICST 

IF(ICST.GE.I)  READ(7)  ICST  .NT  .  N  .KIN ,  T  ,DT .  VOfi .  U .  V,  VORLB 

IF(I.GT.ICST)  THEN 
IF(l.NE.l)  THEN 
WRITE(6.51) 

NT-0 
T-0.0 
ENDIF 
END  IF 
ICST- I 

IF(DTI .NE.0. )  DT-OTI 

C  //////////TIME  STEP  LOOP//////////////////// 

C  .START  COMPUTATIONS  FOR  SUBSEQUENT  TIME  STEPS 
C  FROM  TIME  LEVEL  NS1  TO  NTMAX  (DO  LOOP  1001) 

DO  1001  LOOP-1. NTMAX 
TO-T 
NTO— NT 
DTOOT 
OT-OT+OTINC 
IF  (DT.GT. .08)  OT-.08 
T-T+OT 
NT-NT+1 
DO  1 10  J=1 , N 
DO  1 10  1=1 , IM 
110  VOROLD( I . J )-VOR( I . J ) 

C  DEFINE  AIRFOIL  MOTION 

I F( ICST . EQ . 2)  THEN 

ALP-(T-SIN(7 . 5«T)/7 . 5) «FR.0 . 5 
OMG— (1 . -C0S(7 . 5*T))»FR»0.5 
OMGD— SIN(7.5*T).FR.0,5.7,5 


144 


uViAi.  uai 


*- ,-**/"*V V  V V *V 


o  o  o  o 


r*v\>  •>  ->  ->  -.*  ■ 


END  IF 

I F ( 1 CST .  EQ .  3 )  THEN 

alp»fr*t 

OMG— FR 

I F(ALP .GT .0.6)  THEN 
ALF>-0 . 6 
OMG-O. 

END  IF 
END  IF 

IF( ICST . EQ . 4)  THEN 
ALP— ( 15.-10. *COS( FR»T) ) »PI/1 80 . 

OMG— 10- •FR»SIN(FR«T). PI/180. 

OMGD— 10. •FR*FR»COS(FR*T)*PI/180. 

END  IF 

C. .VELOCITY  BOUNDARY  CONDITION  ON  BODY 

IF(ICST.GT.I)  THEN 
DO  115  K-1 ,KFC 

AS1(K)-AS2(K)*OMG 
BS1 (K)«BS2(K)»0MG 
CS1(K)-CS2(K)»OMG 
DS1(K)-OS2(K).OMG 
115  CONTINUE 
END  IF 

ALPD-ALP»180./PI 

UI-COS(ALP) 

VI~SIN(ALP) 

IF  ( INT(ALPD*10. )  .EQ.  I NT (ALLP(NALLP)»1 0. ) )THEN 
ICPL  -  0 
ICLD  -  0 
ICOUT  -  0 
NALLP  -  NALLP+1 

WRITE(6.61)  NT.T.DT.ALPD.OMG.OMGD 
ELSE 

ICPL  -  1 
ICLD  -  1 
ICOUT  -  1 
END  IF 

ICPL-MOD(NT.NTPLl 
ICLD-MOD(NT.NTLO) 

I COUT-MOD  (  NT  .  NTOUT  ) 

WRITE(6.61)  NT.T.DT.ALPD.OMG.OMGD 

C.. SOLVER  FOR  KINETICS 

CALL  KNTCS(WMIN, URBI .URBP . KMAX , DRMX , KODE , KKK .ALP) 
IF(KOOE.EQ.I)  GO  TO  5000 
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c. .UPDATE  ARRAY  KIN 

DO  120  1-1. IM 
JI»KEN(I)+3 
IF(JI.GT.JR)  JI-JR 
IF(N.LT.JI)  N-JI 
KIN(I)-JI 
120  CONTINUE 

IF(N.LT.NPL)  N-NPL 

C.. SOLVER  FOR  KINEMATICS 

CALL  KMTCS 
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C.. EVALUATE  VORTICITY  LEAVING  BOUNDARY 


DO  125  1-1 , IM 

IF(K£N(I) .GE.NLB)  VOR LB= VOR LB+O T»RDTET»(V( I , NLBP) «0 . 5» (VOR( I . NLB) 
1  +VOR( I , NLBP) )-VSC» (VOR( I , NLBP)-VOR( I , NLB) )/(R2(NL8P)-R2(NLB) ) ) 
125  CONTINUE 


C . . EVALUTE  THE  TOTAL  VOPTICITY  INCLUDING  SOLID  ROTATION 


WSUM1-2 . *OMG*AF+VORLB 
DO  130  J-1.NLB 

ACA-R2(J)*DTET»R1D(J) 

DO  130  1-1 .IM 

WSUM1-#SUM1+VOR( I . J ) *H( I , J ) *ACA 
130  CONTINUE 

WRITE(6 , 57)  KKK.WSUM1 , VOR LB 


C. .CALCULATE  CP 


IF(  ICLD .  EQ. 0  .OR.  ICPL.EQ.0  .OR.  LOOP . EQ . NTMAX )  THEN 
CALL  CPVAL 
DO  135  1-1 .IM 

I F( I CTUR . EO . 1  .AND.  IOTUR( I ) . EQ. 1 )  THEN 
WB( I )«USTAR(I ) *USTAR( I )/VSC 
I F(VOR( I , 1 ) . LT . 0. )  WB( I )— ~WB( I ) 

ELSE 

W8( I )-VOR( 1,1) 

ENDIF 

CONTINUE 

WB(IM1)-WB(1) 

*RITE(4)  NT ,NTPL, T ,RF ,W8 .CP .OMG.OMGD , ALP 
ENDIF 


C.  OUTPUT  AT  EVERY  NTOUT  TIME  STEP 


IF(ICOUT. EQ.0  .OR.  LOOP . EQ . NTMAX )  THEN 
WRITE(6.63)  (VOR(I.I).I-I . IM) 

WRITE! 6,67)  (U( 1.2), 1—1 . IM) 

WRITE(6 , 69)  (KEN( I ) , 1—1 , IM) 

WRITE(21)  W1.W2 

IF( ICTUR.EQ. 1 )  WRITE(6.70)  IOT.IOT8.IM 
ENDIF 


C. . CALCULATE  STREAM  FUNCTION  AND  VORTICITY  AT  VELOCITY 
C  GRID  POINTS  AND  STORE  ON  TAPE  FOR  GENERATING  PLOTS. 
C  INTEGRATION  IS  FIRST  ORDER  TRAPAZOIDAL  RULE. 


I F( ICPL.EQ.0)  THEN 
DO  140  1-1, IM1 
140  POP( I , 1 )— 0. 


C. INTEGRATE  TANGENTIAL  VELOCITIES 


DO  150  1-1 . IM 
DO  149  J-2.NPL 

WOW(I ,J)-.5*(VOR(I , J)+VOR(I , J-1)) 

POP(I ,J)-POP(I ,J-1)-.5»(U(I . J-1)+U(I .J))»R1D(J-1) 
CONTINUE 

WOW( I . 1 )— 2. *VOR( I . 1)-WOW(I,2) 

CONTINUE 
DO  151  J-1 ,NPL 

WOW(IM1 ,J)-WOW(1 . J) 

POP( IM1 ,J)-POP(1 .J) 

CONTINUE 


WRITE(9)  NT , NPL . T , ALPO , Rf , RE , WOW , POP 
END  IF 

1001  CONTINUE 

C  /////////////////END  OF  TIME  STEP//////////////////////// 

C. .ABNORMAL  EXIT.  WRITE  DATA  FROM  LAST  COMPLETED  TIME  STEP 

5000  IF  (KOOE.EQ.1)  THEN 
WRITE(6.52)  KKK 
WRITE(6.53)  NTO 
WRITE(6 , 63)  (VOR(I.I).I-I.IM) 

WRITE(6,67)  (U( 1.2). 1—1 . IM) 

WRIT£(6,69)  (KEN(I).I-I.IM) 

IF(ICTUR.EQ.I)  WRITE(6 , 70)  IOT.IOTB.IM 
WRITE(8)  ICST.NTO.N.KIN, TO, DTO , VOROLD , U, V, VORLB 
STOP 
ENDIF 

C.. NORMAL  EXIT.  WRITE  THE  REQUIRED  DATA  FOR  NEXT  RUN 
WRITE(6,53)  NT 

WRITE(8)  ICST, NT . N, KIN, T.DT,VOR,U,V, VORLB 
WRITEI11)  IDIM.JOIM 
WRITE(II)  1  .  .ALPO,  1000000.  .1  . 

DO  600  1-1 .IDIM 
DO  600  J-1 , JDIM 
Q1  (I , J)-1 
600  CONTINUE 

WRITE(II)  ((QI(I.J).I-I .IDIM) .J-1 .JDIM). 

1  ((Q2(I.J). 1-1. IDIM), J-1.JDIM), 

2  ((03(1. J). 1-1. IDIM). J-1. JDIM), 

3  ((04(1. J). 1-1, IDIM), J-1, JDIM) 

STOP 

51  FORMAT(//20X. •♦AIRFOIL  MOTION  HAS  BEEN  CHANGED,  TIME  RESET-',//) 

52  FORMAT(/10X , 1 -NO  CONVER.  SOLU.  WITHIN' , 15, 3X ITERATIONS* * ,//) 

53  FORMAT (//. 30X . '$$  NEXT  RUN  WITH  NSTART  -',16,'  $$•) 

54  F0RMAT(//3X. ' INPUT  FILE  TQ  ZONST'/) 

57  FORMAT (12X. 'I TKNT  -',I4,9X.'TV0R  -  ’,£13.6,’  VORLB  -\E13.6) 

61  FORMAT(/2X.90(1H-)/1X,  ’  NT  -• , 14. • . 1 , ■  T«\F8.4.'  DT  , F8. 6 

1  AA  — ' , F7.4, '  OMG  — ’ , F6 . 4 ,  OMGD  - ’ . F6 . 4/2X, 90( 1H-) ) 

63  FORMAT  (//I  H  .  ’VORTICITY  - FIRST  RING - ' 

1  ,'CCW  FROM  TRAILING  EDGE0O'/,/(10E13. 6)) 

67  FORMAT  (/I  H  .'TANGENTIAL  VELOCITY  - J-2 - ' 

1  ,'CCW  FROM  TRAILING  EDGEO0 ’/,/( 10E1 3 . 6) ) 

69  F0RMAT(/1H  . 'VORTICITY  EXTENT  00  -»-  KEN  ’ 

1  ,'•••’ ,/,/(1H  ,4(2X .  1013) ) ) 

70  FORMAT (/I H  .'TURBULENT  REGION  00  '.13,'-  1  AND  '.13. ’-’.13) 

END 
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SUBROUTINE  KNTCS(WMIN,URBI , URBP , KMAX , ORMX , KODE , KKK , ALP ) 


KINETICS  OF  THE  PROBLEM 

CALLS  00  WSURFT  WSURF 
EDDYS  FOCFT 
VORTY 

CALLED  BY  00  ZONST 


PARAMETER  f IDII«W0 . JOIM-60) 

PARAMETER  ( IP1-81 , JP1-61 ) 

PARAMETER  (KFC1-41) 

COMMON/ 10/ LOOP . NT , NTMAX , NTOUT 
COMnCN/VTEXT/KIN( IOIM) ,KEN( IDIM) 

COMMON/SGA/GA( IDIM, JDIM) ,GB( IDIM. JDIM) ,GC( IDIM. JDIM) 
COMJON/RHS/GD ( IDIM,  JDIM)  ,DIP1  (  IDIM ,  JDIM)  .DIM1  ( IDIM,  JDIM) 
C0M40N/*FC/AA( JDIM. KFC1 ) ,B8( JDIM, KFC1) 

COMkON/COE/AF,UI , VI , OMG , VSC , NPL , ICTUR, ICST.ICPL 
C0MJ0N/TRIG/CS(KFC1 , IDIM) ,SN(KFC1 . IDIM) 

COMMON/DEITA/DS , DTET , DT 
COt«WON/ZNS/ 1 B 1  ,  IV1  .  IV2,  IB2,  IV1R 

COt*»N/TURl/USTAR(IDIM) ,  EDDY(  IDIM,  JDIM)  .  IOTUR(IDIM) . 

1  YN(KFC1 , JDIM) , IOT, IOTB 

COfcMON/GRD/IM , IM2 . KFC . JR , N 
COMJON /SCALE/H( IDIM, JDIM) 

COMMON/W/VOR(  IDIM,  JDIM)  ,VOROLD( IOIM,  JDIM) 

C0»O«NAEL/U(  IDIM.  JP1  )  ,  V(  IDIM,  JP1 ) 

C0M40N/RGRD/R1D( JP1 ) . R2( JP1 ) 

C0MM0N/DKIN/A2(JDIM) ,A4(JDIM) ,C2(JDIM) ,C4(JDIM) ,D5( JDIM) 
C0MM0N/FM/GAM( IDIM . JOIM) 

KOOE-0 
IKS-1 

IFfALP.EQ.0.)  IKS-IM2+2 

I F (MOD(NT . NTOUT ) . EQ . 0  .OR.  LOOP . EQ . NTMAX )  WRITE(6,10) 

..COMPUTE  FRICTION  VELOCITY  AND  EDDY  VISCOSITY 

IF( ICTUR. EQ.1)  THEN 
CALL  WSURFT 
CALL  EDDYS 
END  IF 

. .CONSTRUCT  DATA  ARRAYS  THAT  ARE  NEEDED  IN  COMPUTATION  OF 
FINITE  DIFFERENCE  FORM  OF  VORTICITY  TRANSPORT  EQUATION 
TIME  TERM  -  FORWARD  DIFFERENCE 

DIFFUSION  TERMS  -  CENTRAL  DIFFERENCES 
CONVECTION  TERMS  -  2ND  UPWIND  DIFFERENCES 

DO  101  I-IKS, IM 
IPP1=»I  +  1 
IM1-I-1 

IF(IPP1 .GT. IM)  IPP1-IPP1-IM 
IF(IMI.LT.t)  IM1-IM+IM1 
JL-KIN(I) 

C.  DETERMINE  IF  REGION  IS  BL  OR  NS 
ICBL-0 


I F ( ( I . GT . I VI . AND . I . LT . IB1 ) . OR . ( I . GT . IB2 . AND . I . LT . I V2) )  1CBL-1 


C.. CONSTRUCT  MATRIX  COEFFICIENTS  FOR  INTERIOR  VORTICITY  SOLUTION 

DO  100  J-2.JL 

T1-H(I,J)*R2(J)/DT 

GA(I.J)-A*(J) 

GB( I , J)— T1— A4( J )-C4( J ) 

GC(I,J)-C4(J) 

GD(I,J)-T1*VOROLD(I,J) 

C.. TURBULENT  FLOW  TERMS  AND  BOUNDARY  LAYER  TERMS 

IF( ICTUR. EQ. 1  .AND.  IOTUR( I ) . EQ. 1 )  THEN 
JPP1-J+1 

1F(JPP1 .GT. JR)  JPP1-JR 
GA( I , J )-GA( I . J )+A4( J ) *EDDY (I , JPP1 )/VSC 
GB(  I .  J  )^B  ( I .  J  )-(  A4(  J  )+C4(  J  )  )  •  EDDY  (I  ,J)  /VSC 
GC( I , J )-GC( I , J )+C4( J ) »EDDY( I . J-1 )/VSC 
END  IF 

IF(ICBL.EO  1)  THEN 
DIP1(I. J)-0. 

DIM1(I.J)-0. 

ELSE 

DIP1 ( I , J)— 05( J) 

DIM1(I.J)-05(J) 

GB(I,J)-GB(I,J)+2.»05(J) 

IF( ICTUR. EQ. 1  .AND.  IOTUR( I ) . EQ. 1 )  THEN 
DIP1  (I , J)— OIP1 f I , J)+05f J)*EDOY( IPP1 , J )/VSC 
D IM1  ( I . J )-OIM1 ( I . J )+D5( J ) *EDDY( IM1 , J ) /VSC 
GB( I .  J)^5B( I . J)+2. *D5( J)»EDDY( I , J)/VSC 
END  IF 
END  IF 

C. .DISCRETIZATION  OF  CONVECTION  TERM  IN  RADIAL  DIRECTION 

VR-V(I,J+1) 

VL-V(I.J) 

A2VR-A2 ( J ) «VR 
C2VL-C2(J)*VL 
IF(VR. LT.0. )  THEN 
GA(I.J)^A(I.J)+A2VR 
ELSE 

G8(I.J)-GB(I,J)+A2VR 

ENOIF 

IF(VL. LT.0. )  THEN 
GB(I,J)-GB(I,J)+C2VL 
ELSE 

GC(I.J)-GC(I,J)+C2VL 
END  IF 

C. .DISCRETIZATION  OF  CONVECTION  TERM  IN  TANGENTIAL  DIRECTION 

UL-0.25«(U(IM1 . J+1 )+U( IM1 ,J)+U(I,J+1)+U(I,J)) 

UR-0. 25* (U( I . J+1)+U(I.J)+U( IPP1 .J+1)+U(IPP1 .J)) 
DTUL-UL/DTET 
DTUR-UR/DTET 
I F(UL .LT.0.)  THEN 
GB(I.J)-G8(I ,J)-OTUL 
ELSE 

DIM1 ( I , J)— OIM1 ( I , J)+OTUL 
ENOIF 

IF(UR .LT.0.)  THEN 
DIP1 ( I . J)— OIP1 ( I , J)-OTUR 
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G8(I , J)-GB(I . J)+DTUR 
END  IF 

100  CONTINUE 

C. .NEUMANN  TYPE  B.C 

IF(JL.GE. JR)  G8(I ,JL)-G8(I , JL)+GA( I , JL) 

101  CONTINUE 
KKK-0 
ICMR-2 
DMAXP-100. 

UR8-UR8I 

C..VORTICITY  CONVERGENCE  LOOP 

500  CONTINUE 
KKK-KKK+1 

C. . BOUNDARY  CONDITION  FOR  KINETICS  IN  TURBULENCE 

IF( ICTUR . EU. 2)  THEN 
DO  110  1-1. IM 
I F( IOTUR( I ) . EQ. 1 )  THEN 
IN-I 

IF(IN.GT.KFC)  IN-IM-I+2 
YP2-YN  ( I N .  2  )  «UST AR  ( I )  ASC 
I F(YP2 .GT . 5. )  THEN 
IOTUR( I )— 0 
ELSE 

VOR(I.2)-VOR(I,1) 

END  IF 
ENDIF 

110  CONTINUE 
ENDIF 

C.. SOLVE  FOR  EACH  ZONE;  USE  ONLY  HALF  OF  GRID  IF  SYMMETRIC  FLOW 

IF(ALP.NE.0.)  THEN 
CALL  VORTY(IB1 . IB2, 1 , ICS.0.KOOE) 

I F ( ( IV1+1 ) . LT . IB1 )  CALL  VORTY( IB1-1 , IV1 ,-1 . ICB1 , 1 .KODE) 
IF( IV2. GT. IB2)  CALL  VORTY( 182+1 , IV2 , 1 , ICB2 . 1 .KODE) 

CALL  VORTY( IV2.IV1R.1 ,ICV, 0.KOOE) 

ELSE 

CALL  VORTY(IM2+2.IB2.1, ICS. 0, KODE) 

CALL  VORTY(IB2+1 .IV2.1  ICB2.1 .KODE) 

CALL  VORTY(IV2, IM, 1 . ICV.0.KODE) 

DO  120  1-2. IM2 
II-IM+2-I 
DO  120  J-2.N 

120  VOR(  I .  J )— VOR(  1 1  ,  J  ) 

EM)I  F 

IF(KODE.GT . 0)  RETURN 

C.. DETERMINE  VOR.  EXTENT  AT  EACH  RADIAL  LINE  IN  N-S  REGION 
C  AND  STORE  IN  ARRAY  KEN 

DO  130  1-1 . IM 
DO  129  J— N ,  2 .— 1 

129  IF(A8S(VOR(I  .-0)  GT.WMIN)  GO  TO  130 

130  KEN( I )— J 

C.. DETERMINE  FOURIER  COEFFICIENTS  OF  VORTICITY  AND  EVALUATE 
C  SURFACE  VORTICITY  BY  UNDER-RELAXATION  TECHNIQUE 


DO  140  1-1 . IM 
JL-KEN(I) 

JLP-JL+1 
DO  139  J-1 . JL 

139  GAM( I , J )«H( I , J )*VOR( I , J ) /FLOAT ( IM2) 

IF(JLP.LE.N)  THEN 

DO  138  J-JLP.N 
1 38  GAM( I , J)-0. 

END  IF 

140  CONTINUE 
CALL  FOCFT 
CALL  WSURF 

MI-0 
MW  1-0 
INMAX— 0. 

DMAX— 0. 

W1-0. 

DO  150  1-1, IM 
W1— AA( 1 , 1 )•  .  5 
DO  149  K-2.IM2 

149  W1-W1+AA(1 ,K)*CS(K,I)+BB(1  ,K)*SN(K,I) 

W1-**1+AA(1  , KFC)*CS(KFC , I )•  .  5 
W1«W1/H( 1,1) 

UR1-URB»H( I , 1 )«»UR8P 
1WMW1  *UR1-f(  1  .-UR1  )*VOR(  1,1) 

AVWWkBS(WW) 

IF(AWW.LT.10.)  AWW-10. 

DW-ABS (WW-VOR (1,1)) /AWW 
IF(DW.GE.DMAX)  OMAX-OW 
IF(DW.GE.DMAX)  MI-I 
IF(AWW.GE.WMAX)  WMAX-AWW 
IF(AWW.GE.MMAX)  MWI-I 
VOR(I.1)-WW 

150  CONTINUE 

C. .ADJUST  UNDER-RELAXATION  PARAMETER  URB. 

C  EXIT  KNTCS  IF  CONVERGED.  CONTINUE  ITERATIONS  IF  NOT. 

C  ABORT  IF  MAXIMUM  ITERATIONS  EXCEEDED. 

IF(KKK.EO. ICMR)  THEN 
DDM- ( DMAX P— DMAX ) /DMAX 
IF(DDM. LT .0 . )  URB-0 . 90»URB 
IF(DDM.LT.0. 1 .ANO.ODM.GT.0. )  URB-1 ,08*URB 
IF(DDM.GT.0.1)  URB-1. 05*URB 
IF(URB -GT .0 . 65)  URB-0. 65 
ICMR-ICMR+2 
END  IF 

DMAXP—DMAX 

IF(MOO(NT,NTOUT) -EO.0  .OR.  LOOP . EQ.NTMAX) 

+  WRITE(6, 12)KKK.DMAX,MI , VOR(MI . 1 ) .WMAX.WNI ,URB, ICS, ICB1 , ICB2 , ICV 
IF(DMAX. LE.DRMX)  RETURN 
IF(KKK. LT.KMAX)  GO  TO  500 

C.  NO  CONVERGENCE  OCCUR ED;  RETURN  TO  MAIN  WITH  KODE-1 
KOOE-1 

IF(MOO(NT,NTOUT) .NE.0  .OR.  LOOP .NE. NTMAX)  THEN 
WRITE(6, 10) 

WRITE(6,12)KKK,0MAX,MI  ,VOR(MI  ,  1  )  .VWAAX  .MMI  ,  URB  ,  ICS,  1CB1  ,  ICB2.ICV 
END  IF 
RETURN 

10  FORMAT (//.25X, 'SURFACE  VORTICITY  INFORMAT IONOO ’ . 20X , 

1  'INTERIOR  VORTICITY  INFORMATIONS '/'  ITER',6X. 
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SUBROUTINE  WSURFT 


EVALUATE  BOUNDARY  VORTICITIES  IN  TURBULENT  FLOW  REGION 


CALLS  OO  USTARF 
CALLED  BY  OO  KNTCS 


PARAMETER  (IDII*W0. JDIM-60) 

PARAMETER  (KFC1-41) 

COMMON/SCALE/H( IDIM, JDIM) 

COMMON/W/VOR ( IDIM.JDIM)  ,  VOROLD(  IDIM.JDIM) 

COMMON/TUR1/USTAR( IDIM) , EDDY( IDIM, JDIM) , IOTUR( IDIM) , 

1  YN(KFC1 . JOIM) , IOT , IOTB 

COhMON/COE/AF .UI,VI,OMG,VSC.NPL.ICTUR,ICST,ICPL 
COMMON/GRD/IM . IM2 . KFC . JR , N 

.EVALUATE  USTAR 

OMG2-2 . »OMG 
USTAR(1)-0. 

DO  100  1-2. IM 
IN— I 

IF(IN.GT.KFC)  IN-IM-I+2 
G»*»  (  VOR  ( 1 . 1  )-OMG2  )  •  YN  ( I N .  1 ) 

AGWM^BS(GM) 

YP-SQRT(AGM*YN(IN,1)/VSC) 

US7A.T(  I )— AGM/YP 

IF(YP . GT . 5. )  USTAR( I )-USTARF(USTAR( I ) . AGM . YN( I N . 1 ) . 5 . . -3 . 05 . YP) 
IF(YP.GT.3B.)  USTAR(I)-USTARF(USTAR(I).AGM,YN(IN,1).25.5.5.YP) 
100  CONTINUE 
RETURN 
END 

FUNCTION  USTARF (USTR . UV , YN , A . C , YP) 


EVALUATE  USTAR  ITERATIVELY  IN  INERTIAL  LAYER 
CALLS  00  NONE 
CALLED  BY  00  WSURFT 


COMMON/COE/AF , U I .VI ,OMG, VSC.NPL. ICTUR. ICST. ICPL 
DO  100  1-1 .20 
YP— YN»USTR/VSC 
USTARF-UV/ ( A  » A  LOG ( YP ) +C ) 

ABSER— ABS((USTARF-USTR)/USTARF) 

I F( ABSER.LT. 0.005)  RETURN 
USTR-USTARF 
100  CONTINUE 

WRITE(6. 10)  YP.  USTR, ABSER 

10  F0RMAT(2X , 'NO  CONVERGENCE  IN  USTARF.  YP , USTAR , ABSEROO 1 , 3F9 . 5) 
STOP 
END 
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SUBROUTINE  EDDYS 


EVALUATE  EDDY  VISCOSITY 
CALLS  00  NONE 
CALLED  BY  00  KNTCS 


PARAMETER  ( IDII*W30 ,  JD  11^60) 

PARAMETER  ( IP1-81 . JP1-61 ) 

PARAMETER  (KFC1-41 ) 

PARAMETER  (INOR1-20) 

COMMON/SCALE/H( IDIM, JDIM) 

COMMON/ 10/ LOOP , NT . NTMAX , NTOUT 
COfcMON/W/VOR(IDIM,  JDIM)  ,VOROLD( IDIM,  JDIM) 
C0*iH0N/TUR1/USTAR(  IDIM)  .EDDY(  IDIM.  JDIM)  ,  IOTUR(  IDIM)  . 
1  YN(KFC1 . JDIM). IOT. IOTB 

COfcMON/TUR2/ I NOR ( I NOR 1 .JDIM) . IWK( JDIM. JDIM) 
COAWON/GRD/IM,  IM2 . KFC ,  JR ,  N 

C0IA10N/C0E/AF .UI.VI.OMG.VSC.NPL.ICTUR.ICST.ICPL 
COM*ON/VTEXT/KIN(IDIM) .KEN(IDIM) 

COMHON/VEL/U( IDIM. JPt ) . V( IDIM, JP1 ) 

UC-IM+2 
DO  105  1 1—2. IM 
I-II 

IOTUR( I )-0 
FMAX-O. 

YMAX-0. 

VDMX-0. 

JL-KIN(I)+1 
IFfJL.GT. JR)  JL-JR 
IF(KIN( 1+1 )  .GT .  JL)  JL-KINII+1) 

IF(KIN(I-1).GT.JL)  JL-KIN(I-l) 

C.. CALCULATE  TERMS  INDEPENDENT  OF  Y  AND  INNER  VISC. 

C  INCLUDE  CORRECTION  FOR  THE  NORMAL  DIRECTION 

DO  100  J-2.JL 

I F ( I I . LT . IN0R1+1 )  I-INOR(II.J) 

IF(II.GT.JPI)  I-IMB-INOR(IMB-II.J) 

IN— I 

IF(I.GT.KFC)  II+-IMB-I 
AVOR»ABS ( VOR ( I . J) ) 

YP-USTAR( 1 1 ) *YN( IN . J)/VSC 
IF(YP .GT . 500 . )  THEN 
TI-YN(IN.J) 

ELSE 

T1«YN(IN, J)»(1 ,-EXP(-YP/26. )) 

END  IF 

EDDY(I,J)«0.16»T1*T1*AVOR 
FY=AV0R*T1 
IF(FY.GT.FMAX)  THEN 
FMAX-FY 
YMAX-YN(IN, J) 

END  IF 

UV-0.5*(U(I,J+1)+un.J)) 

W— 0 . 5*(V(  I  .  J+1  )+V(  I  ,  J  )  ) 

VSPHY»(UV*UV+W*W)/H(  I  ,J) 

I F ( VSPHY . GT . VDMX )  VDMX-VSPHY 
100  CONTINUE 


C.. OUTER  EDDY  VISCOSITY 


IF(FMAX . £Q. 0 . )  GO  TO  105 

FWAKE-FMAX.YMAX 

CFWAKE-0 . 25*YMAX*VDMX/FMAX 

IF(CFWAKE.  LT.FWAKE)  FWAKE-CFWAKE 

I C ED-0 

DO  103  J-2.JL 

IF(  I 1 . LT . INOR1+1 )  I-INOR(II.J) 

XF( 1 1 .GT. JP1)  I-IMB-INOR(IMB-II.J) 

IN- 1 

IF(I.GT.KFC)  IhMMB-I 

FKLEB-1 ./(I .+5 .5*(0 . 3«YN( IN , J )/YMAX)»*6 . ) 
CEDDY-0 . 01 68«1 . 1 *FWAKE»FKLEB 
IF(ICED. EQ.0  .ANO.  CEDDY. LT . EDDY( I , J) )  THEN 
ICED-1 

IF(CEDDY/VSC.GT0.)  IOTUR(II)-1 
END  IF 

IF(ICED.EQ.I)  EDDY(I,J)^EDDY 
103  CONTINUE 
105  CONTINUE 

C . . MAKE  SURE  FLOW  REMAINS  TURBULENT  AFTER  TRANSITION 

IOT-0 

TOTB-0 

DO  110  I— KFC, 1 1 

TITmTUR-  t 

IFfIB.GT. IM)  IB-IM 

IF(I0TUR( I ) . EQ. 1  .AND.  IOT.EO.0)  IOT-I 
IF(IOT.NE.O)  IOTUR( I )-1 

IF(IOTUR( IB) . EQ. 1  .AND.  IOTB.EO.0)  IOTB-IB 
IF(IOTB.NE.0)  I0TUR(IB)-1 
110  CONTINUE 

C.  ASSIGN  EDDY  VISCOSITY  IN  THE  WAKE  ALONG  WAKE  GRID 

DO  120  J»^2. JR-1 
I»MN0R(2.JW) 

IBKNIMB— IW 
DO  119  J-JW.JR 
I— IWK( JW. J) 

IF( I , EQ. IN0R(2 . J ) )  GO  TO  116 

115  IB-IMB-I 
IF(I.EQ.I)  IB-1 
EDDY( I , J )— EDOY( IW, JW) 

EDDY( IB , J )— EDOY( I8W . JW) 

116  1—1—1 

IF(IWK(JW-1 ,J) .IT. I)  GO  TO  115 

119  CONTINUE 

120  CONTINUE 
RETURN 
END 


u 


SUBROUTINE  VORTY( IS , I L , INC, IC, ICBL.KODE) 

C 

C  •  *« . 

C  CALCULATE  VORTICITY  BY  USING  LINE- 
C  RELAXATION  METHOD  ON  EACH  RADIAL  LINE 
C 

C  CALLS  00  TRID 
C 

C  CALLED  BY  00  KNTCS 

C  . . . . . . 

c 

PARAMETER  ( IDIVN80 , JDIK*60) 

PARAMETER  (KFC1-41 ) 

COMHON/COE/AF.UI , VI ,OMG ,VSC,NPL , ICTUR , ICST, ICPL 
C0MM0N/TUR1 /USTAR( IOIM) , EDDY( IDIM, JDIM) . IOTUR( IDIM) , 

1  YN(KFC1 . JDIM) , IOT , IOTB 

COfcMON/GRD/IM, I M2 , KFC , JR , N 
COMMON/VTEXT/KIN( IDIM) ,K£N( IDIM) 

COfcMON/UNF/URR , URF , DFMX , NCC 

COfcMON/SGA/GA( IDIM. JDIM) ,GB( IDIM, JDIM) ,GC( IDIM, JDIM) 
COMMON/RHS/GD( IDIM, JDIM) ,DIP1 ( IDIM, JDIM) ,DIM1 ( IDIM, JDIM) 
COMy(ON/VV/VOR(  IDIM.  JDIM)  ,VOROLD(  IDIM,  JDIM) 
COMMON/SOLV/SUB (JDIM) ,OIAG( JDIM) , SUP( JDIM) , RHS( JDIM) 

IF(KOOE.EQ.I)  RETURN 
DO  100  IC— 1  ,  NCC 
VYMAX-0. 

DMAX-0, 

DO  110  II-IS.IL.INC 
I— 1 1 

IPP1-I+1 

IM1-I-1 

IF(I.GT.IM)  I-II-IM 
IF(IPPI.GT.IM)  IPP1-IPP1-IM 
I F( IM1 ,GT . IM)  IM1-IM1-IM 
JI-2 

I F( ICTUR . EQ. 2  .AND.  IOTUR(I)  E0.1)  JI-3 

JIM-JI-1 

JL-KIN(I) 


! 

! 
M  * 


V. 


y 

/ 
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C.. CONSTRUCT  TRIDIAGONAL  MATRIX  AND  SOLVE 

DO  120  J-JI.JL 
J1-J-JIM 
SUP(J1)-GA(I,J) 

DIAG(J1)-GB(I ,J) 

SUB(J1)-GC(I. J) 

RHS(J1)-GD(I , J)+DIP1(I , J)«V0R(IPP1 , J)+DIM1(I , J)«V0R(IM1 . J) 
120  CONTINUE 

RHS(1)«RHS(1)-SUB(1)«V0R(I .JIM) 

CALL  TRIO(JI) 

C . . UNDER-RELAX  THE  RESULT  OF  THE  MATRIX  SOLUTION 

DO  130  J-JI.JL 

IF( ICBL  .EQ.  1)  THEN 
W#-RHS(J-JIM) 

ELSE 

«WMJRF»VOR( I . J)+URR«RHS(J-JIM) 

WABS-ABS(WW) 

DD-ABS(WW-VOR( I , J)) 

IFfDD.GE.DMAX)  DMAX-DD 
IF(WABS.GT.WMAX)  WMAX-WABS 
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k^i  wk-Jlk  ^  Sh  1-. 


V.! 

■5 


V? 


VOR(I ,J)«WW 
130  CONTINUE 
110  CONTINUE 

C..IF  BOUNDARY  LAYER  ZONE,  RETURN  TO  KNTCS  AFTER  SOLVING  EXPLICITLY 
C  ELSE  ITERATE  FOR  CONVERGENCE  AND  RETURN  WHEN  CONVERGED  OR  WHEN 
C  MAX  NUMBER  OF  ITERATIONS  IS  EXCEEDED. 

IF(ICBL.EQ.I)  RETURN 
IF(DMAX/YfMAX .  LE.DFMX)  RETURN 
100  CONTINUE 

WRITE(6, 10)  II.DMAX 

KOOE-1 

RETURN 

10  FORMAT(2X. ’NO  TRIO  CONVERGENCE  AT  I-'. 14,-  DMAX-’ . E10. 3) 

END 


PARAMETER  (JDIM-60) 

COMy©N/SOLV/SUB( JDIM) ,DIAG(JOIM) .SUP(JDIM) .RHS(JDIM) 

DO  100  J-2.N 

RAT  10— SUB  (  J  )/D  I  AG(  J- 1 ) 

DIAG(J)-0IAG(J)+RAT10*SUP(J-1) 

RHS( J )-RHS( J )+RAT IO*RHS( J-1 ) 

100  CONTINUE 

RHS(N)«RHS(N)/DIAG(N) 

DO  110  J-N-1 .1 ,-1 

110  RHS(J)«(RHS(J)-SUP(J)»RHS(J+1))/DIAG(J) 

RETURN 

END 
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SUBROUTINE  FOCFT 


DETERMINE  FOURIER  COEFFICIENTS  AT  EACH  RING 

CALLS  00  NONE 

CALLED  BY  00  VORTY 
CPVAL 


PARAMETER  ( IDII*-80.  JDII**60) 

PARAMETER  (KFC1-41) 

COMk»N/WFC/AA( J0IM.KFC1 ) ,BB(JD1M.KFC1 ) 
C0M*0N/SMT/SGMA(KFC1  ) 

COMrfON/SCALE/H(  ID1M,  JDIM) 

COMMON/GRD/ I M , IM2.KFC. JR.N 
COMMON/VTEXT/KIN( IDIM) ,KEN( IOIM) 
C0M*DN/TRIG/CS(KFC1 . IDIM) .SN(KFC1 . IDIM) 
COfcMON/FM/GAM(  IDIM,  JDIM) 

DO  100  J-1.N 
DO  100  K«1 ,KFC 
BB( J ,  K)«0 . 0 
AA(J,K)«O.0 
100  CONTINUE 

DO  110  1-1. IM 
NI-KEN(l) 

DO  109  J-1 , NI 

AA( J . 1 )— AA( J . 1 )+GAM( I .J) 

DO  109  K-2.KFC 

AA(J.K)-AA(J.K)40AM(I.J)*CS(K.I) 
B8(J.K)-eB(J.K)4GAM(I  ..O.SNfK.I) 

109  CONTINUE 

110  CONTINUE 

00  120  J-2.N 
DO  120  K-2.KFC 

AA(J.K)-AA(J.K).SGMA(K) 

BB(J ,K)— BB( J ,K)«SGMA(K) 

120  CONTINUE 
RETURN 


SUBROUTINE  WSURF 


PARAMETER  (IDIM-80,  JDIHA-60) 

PARAMETER  (JP1-61) 

PARAMETER  (KFC1-41 .KFC2-42) 

COMMON/WFC/AA( JDIM.KFC1 ) ,BB( JDIM.KFC1 ) 

C0M40N/SMT /SGMA ( K  FC 1 ) 

COMn(ON/COE/AF,UI  ,VI  .OMG.VSC.NPL.  ICTUR,  ICST,  ICPL 
COMMON/ ABCOS/ AS 1 (KFC1 ) . BS1 (KFC1 ) ,CS1 (KFC1 ) , DS1 (KFC1 ) 
COMMON/GRO/IM. IM2.KFC, JR.N 
COMMON/COR/RP ( J  P 1 ,KFC2) ,RL(JP1 ) 

COfcMON/VLB/VORLB , NLB 

COfcMON  01(IDIM,JDIM),Q2(IDIM.JDIM)>03(IDIM,JDIM),04(IDIM,JDIM) 

C.. DETERMINE  FOURIER  COEFFS.  OF  BOUNDARY  VORTICITY.  AA( 1 ,M) , BB( 1 ,M) 

C  FOR  M.GE.4 

PI-4.*ATAN(1 . ) 

DO  100  K-4.KFC 
CW1-0. 

DW1-0. 

DO  101  J-2.N 

W1-1 ,/RP(J ,K— 3)— 1 ,/RP(J+1 ,K-3) 

CW1-CW1+AA(J.K)**1 
DW1-0*M+BB(J.K)«W1 
101  CONTINUE 

W2-1 ./(I ./RPM . K— 3)— 1 ,/RP(2 . K— 3) ) 

BB( 1 .K)“W2*( (AS1 (K)— OS1 (K) )«FLOAT(K-3)-OW1 ) 

AA( 1 ,K)— W2*((CS1 (K)+BS1 (K) ) «FLOAT(K-3)+CW1 ) 

100  CONTINUE 

C. .DETERMINE  AA( I.K).BB(I.K),  FOR  K.LE.3 

CW-0. 

CW1-0. 

CW2-0. 

DW1-0. 

DW2-0 . 

DO  130  J-2.N 

CW1-CW1+AA(J ,2)«(RP(J+1 , 1 )-RP( J ,  1  )) 

DW1-DW1+BB(J.2)»(RP(J+1 . I)-RP(J.I)) 

CW2-CW2+AA f J . 3 ) • ( R L ( J+ 1 ) -R L ( J ) ) 
DVK2«DW2+BB(J,3)«(RL(J+1)-RL(J)) 

130  CONTINUE 

DO  140  J-2.NLB 

140  CW"CW+AA( J , 1 ) • (RP( J+1 , 2)— RP( J . 2) ) 

AA(1 ,1)— (4..AF/PI*0MG+CW+2.«V0RLB/PI)/(RP(2,2)-RP(1 ,2)) 

AA(1 . 2)"(2 . *VI— CS1 ( 2 )  — BS 1 (2)-CW1 )/(RP( 2 , 1 )-RP( 1 , 1 )) 
BB(1.2)«(AS1(2)-DS1(2)-2.»UI-DW1)/(RP(2, 1)-RP(1 , 1)) 

AA(1 .3)— (CS1(3)+BS1(3)+CW2)/RL(2) 

BB( 1 .3)-(AS1 (3)-OS1 (3)-OW2)/Rl(2) 

DO  150  K-2.KFC 


AA(1 ,K)«AAf1 , K ) • SGMA (  K  ) 
BB( 1 ,K)-8B(1 . K ) • SGMA (  K ) 


150  CONTINUE 


v  ■--.vS'Cs*: 


ooooooooo 
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SUBROUTINE  KMTCS 

KINEMATICS  OF  THE  PROBLEM 
CALLS  DO  NONE 
CALLED  BY  OO  ZONST 


PARAMETER  ( ID  11^80 ,  JDIV*60) 

PARAMETER  ( IP1-81 . JP1-61 ) 

PARAMETER  (KFC1-41 .KFC2-42) 

DIMENSION  VTOTPH(IDIM.JDIM) . VTOTCO( IDIM. JDIM) 

COMMON/COE/ AF , U I . V I , OMG . VSC , NPL, ICTUR, ICST , ICPL 
COMMON/VTEXT/KIN(IDIM) .KEN(IDIM) 

COMMON/SMT/SGMA(KFC1 ) 

COMMON/ ABCDS/ AS 1 (KFC1 ) ,BS1 (KFC1 ) .CS1 (KFC1 ) , DS1 (KFC1 ) 
COMMON/RHS/A(JP1 ,KFC1),B(JP1 , KFC1 ),C(JP1,KFC1),D(JP1 . KFC1 ) 
COMyCN/*FC/AA(JDIM,KFC1  )  ,BB(  J0IM.KFC1  ) 

COMMON/COR/RP( JPt  ,KFC2)  ,RL(JP1) 

C0M40N/VEH/UC( IP1 , JP1  ) , VC( IP1  , JP1  ) 

COMMON/ W/VOR(  IDIM.  JDIM)  ,  VOROLD(  ID  IM,  JDIM) 

COMKONAEL/U(  IDIM.  JP1  )  .  V(IDIM, JP1 1  ,HSTAR(  IDIM.  JDIM) 

COMMON/GRD/ 1 M .  I  M2  .  KFC  ,  JR  .  N 

COMMON/TR I G/CS ( K  FC 1 , IOIM) .SN(KFC1 . IDIM) 

COMMON  01 (IDIM. JDIM). Q2( IDIM, JDIM), Q3( IDIM, JDIM), Q4( IDIM. JDIM) 
COMPLEX  HSTAR , VTOTCO . VTOTPH 

NP-N4-1 

DO  100  K-1 , KFC 


A< 

1.K 

-AS1| 

K) 

B< 

1.K 

-esii 

K 

c( 

1.K 

-CS1< 

K 

D( 

1.K 

-0S1( 

K) 

100  CONTINUE 

DO  500  J-2.NP 
JM1-J-1 


C.. BOUNDARY  VELOCITY  CONTRIBUTIONS 
W1-.5/RP(J.2) 

C(J,2)«W1»(C(1 .2)— B( 1 . 2 ) ) +V I 
D(J.2)-N1*(D(1 .2)+A(1 ,2))-UI 
A(J,1)-A(1,1)/RP(J.1) 

A(J,2)-W1»(0(1 ,2)+A(1 .2) )+UI 
B( J ,2)»W1 *(B( 1 ,2)-C( 1 , 2) )+VI 
DO  110  K-J.KFC 
W1-.5/RP(J,K) 

C(J,k)-M1»(C(1 .K)-B(1 ,  K)  ) 

D(J.K)-N1»(0(1 . K ) +A ( 1  . K )  ) 

A(J ,K)«D(J ,K) 

B(J.K)— C(J.K) 

1 10  CONTINUE 

C. . INTEGRATE  ALONG  RADIAL  AXIS  FOR  AA(J.K),  BB(J,K).  (K.GE  4) 


DO  120  K-4.KFC 


OW1-0. 

00  119  JJ-1 . JM1 

WI^PlJJ+l  ,K+1  )— RP(  JJ  ,K+1  ) 

CW-AA(JJ,K)»W1+CW 

D#«eB(JJ,K)«W1+DW 

119  CONTINUE 
IF(J.LT.NP)  THEN 

DO  118  JJ-J.N 

W1-1 ,/RP(JJ.K-3)-1 ./RP(JJ+1 , K-3) 
CV¥1-AA(JJ,K)*W1+CW1 
D#1-8B(JJ,K)*W1+DW1 
118  CONTINUE 
ENDIF 

W2- . 5/( FLOAT ( K+1 ) «RP ( J , K ) ) 

CW2«*2*CW 

0W2HW2*0W 

W3-.5»RP(J,K-2)/FL0AT(K-3) 

CW3-W3*CW1 

0W3«W3*DW1 

C(J.K)-C(J,K)+CW2-CW3 
D( J ,K)-0( J ,K)+0W2-0N3 
A(J ,K)« A( J ,K)+0W2+0W3 
B( J  .K)-B(J .K)-CW2-CW3 

120  CONTINUE 

C. .COMPUTE  AA(J,K),  BB(J.K);  (K.LE.3) 


CW-0. 

CW1-0. 

CW2-0. 

CW3-0. 

CW4-0. 

DW1-0. 

DW2-0. 

OWJ-0. 

DW4-0. 

DO  130  JJ-1 , JM1 

CV¥-C»+AA(JJ,1)*(RP(JJ+1  . 2)-RP(  J  J  ,  2)  ) 
CW1«C01+AAlJ  J ,2)«(RP(JJ+1 , 3)-RP( J J , 3) ) 
DW1^X»1+BB(JJ.2)*(RP(JJ+1  . 3)-RP(  J  J ,  3 ) ) 
CW2^N2+AAl  J  J  ,3)*(RP(  JJ+1  , 4)-RPf  J  J  ,  4) ) 
DW2HJW2+88(  JJ ,3)*(RPl JJ+1 , 4)-RP( JJ , 4) ) 
130  CONTINUE 

IF(J.LT.NP)  THEN 
DO  140  JJ-J.N 

CW3-CW3+AA( J J . 2 )•( RP ( J J+ 1 . 1 )-RP ( J J . 1 ) ) 
DW3— 0W3+BB ( J  J ,2)*(RP(J J+1 . 1 )-RP(JJ , 1 )) 
CW4-CW4+ AA ( J J . 3 } • ( R L ( J J + 1 ) -R L ( J J ) ) 

D#4— 0W4+8B(  J  J  ,3)*(RL(JJ+1)-RL(JJ)) 

140  CONTINUE 
ENDIF 

C( J . 1 )— (C( 1.1)+. 5*CW)/RP( J . 1 ) 

CW3— .  5*CY13 
DW3-.5»0W3 
W2-1 ,/(6..RP(J,2)) 

CW1-CW1»#2 
DN1-0W1  »\*2 

C(J  ,2)><(J  ,2)+CVM-C*3 
D(J.2)-0(J,2)+0W1-0#3 
A(J.2)-A(J,2)+0#140W3 
8(J,2)*43(J  .  2)-CW1-CW3 
#2-1 ,/(8.*RP(J.3)) 

CV*2’<W2»W2 

DW2«0W2*W2 
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W2-. 5»RP( J  .  1 ) 

CW4-CW4*W2 

0W4-DW4.W2 

C(J.3)-C(J.3)+CW2-CW4 
D ( J . 3 ) -0 ( J . 3 ) +0W2-0W4 
A(J,3)-A(J.3)+DW2+DW4 
B(J ,3)-8(J ,3)-CW2-CW4 
500  CONTINUE 

DO  150  K-2.KFC 
DO  150  J-2.N 

A(J.K)-A(J,K)*SGMA(K) 

B(J,K)«8(J , K)»SGMA(K) 

C(J.K)-C(J,K)*SGMA(K) 

D(J,K)«0(J,K)«SGMA(K) 

150  CONTINUE 

C.. ABOVE  COEFFICIENTS  ARE  FOR  V(RHO)  *  V(PH1)  IN  NON-ROTATING 
C  FRAME  OF  REFERENCE  IN  COMPUTATIONAL  PLANE. 

C  CALCULATE  VELOCITIES  WITH  THESE  VALUES. 

DO  160  1-1 .IM 
JL— KIN(I  )+1 

I F ( ICPL.  EQ. 0  .AND.  NPL.GT.JL)  JL-NPL 

IPP1-I+1 

IM1-I-1 

IF(I.EQ.I)  IM1-IM 

IF(I.EQ.IM)  IPP1-1 

IF(KIN( IPP1 ) .GE. JL)  JL— KIN( IPP1 )+1 

IF(KIN(IM1) .GE.JL)  JL-KIN( IM1 )+1 

DO  159  J-2.JL 

U(I.J)-C(J.1)..5 

V( I ,  J )— A( J , 1 ) * . 5 

DO  158  K-2.IM2 

U(I,J)-U(I. J)+C(J,K).CS(K.I)+0(J,K).SN(K.I) 
V(I.J)-V(I1J)+A(J.K).CS(K.I)+B(JIK).SN(K,I) 

158  CONTINUE 

U(I.J)— U(I,J )+C( J ,KFC) «CS(KFC .  I )•  . 5 
V(I.J)-V(I,J)+A(J.KFC).CS(KFC,I)*5 

159  CONTINUE 


C.. TRANSFORM  VELOCITY  COMPONENTS  TO  ROTATING  FRAME  OF  REFERENCE 
C  FOR  USE  IN  KNTCS 

DO  157  J-2.JL 

W1»V(I , J) «CS(2, I)-U(I,J)»SN(2,I) 

W2-V(I , J)*SN(2,I)+U(I ,J).CS(2. I) 

W1-W1+0MG«UC( I .J) 

W2-W2+0MG*VC(I , J) 

IF(  J.LE.60  )THEN 
VTOTCO(I.J)  =  CMP LX  (W1 ,W2) 

VTOTPH(I.J)  -  VTOTCO( I . J )/HSTAR( I . J ) 

02(1. J)  -  REAL(VTOTPH( I . J ) ) 

Q3( I . J )  -  AIMAG(VTOTPH( I , J)) 

END  IF 

U( I , J)— W2»CS(2, I)-W1*SN(2, I) 

V( I , J)-W1 «CS(2 . I )+W2«SN(2 , I ) 


157  CONTINUE 
160  CONTINUE 


oooooooooo 
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SUBROUTINE  CPVAL 


CALCULATE  CP  VALUES  BY  CP  INTEGRAL 

CALLS  00  COED 
FOCFT 

CALLED  BY  00  ZONST 


PARAMETER  ( ID1M-80 , JDIM-60) 

PARAMETER  ( IP1-81 . JP1-61 ) 

PARAMETER  (KFC1-41) 

PARAMETER  (1 CP-40) 

COMMON/GRD/IM, IM2.KFC, JR,N 
C0MM0N/TRIG/CS(KFC1 , IOIM) , SN(KFC1 . IDIM) 

COMMON/ WAOR(  IDIM,  JDIM)  .VOROLD(IDIM.  JDIM) 

COMMON/COE/AF , U I ,VI ,OMG,VSC,NPL. ICTUR , ICST . ICPL 
C0MM0N/SMT/SGMA(KFC1 ) 

C0M40NAEL/U(  IDIM, JP1  )  ,V(IDIM,  JP1  ) 

COMMON/VTEXT/KIN( IOIM) .KEN(IOIM) 

C0MM0N/CPC/CCP(KFC1 , JDIM) ,CP( IP1 ) 

C0Ml40N/WFC/AA(JDIM.KFC1 ) ,BB(JDIM,KFC1 ) 

COMMON/RHS/ A ( J P  t .KFC1) ,B(JP1 ,KFC1) ,C(JP1 ,KFC1) ,D(JP1 , KFCt ) 
COMMON/ FM/GAM ( ID IM , JDIM) 

COMMON/CPW/AV*(  IDIM)  ,BW(  IDIM)  ,CW(  I  CP)  ,D#(  ICP)  .  EW(  ICP)  .  FW(  ICP)  , 
t  AA2 ( I CP ), BB2 ( I CP ). CC2 ( I CP ). DD2 ( I CP ) 

COMMON  01 (IDIM, JDIM) ,Q2( IDIM, JDIM) ,Q3( IDIM, JDIM) ,04( IDIM. JDIM) 
DIMENSION  UCC( JDIM.KFC1 ) ,UCS(JDIM,KFC1 ) . VCC( JDIM.KFC1 ) , 

1  VCS( JDIM.KFC1 ) , FC(KFC1 ) . FS(KFC1 ) 

DO  105  1-1, IM 
JL-KEN(I) 

JLP-JL+1 
DO  100  J-1 , JL 

100  GAM( I , J)— VOR( I , J)/FLOAT ( IM2) 

IF(JLP.LE.N)  THEN 
DO  103  J— JLP.N 
103  GAM( I . J )— 0 . 

ENDIF 

105  CONTINUE 

DO  110  1-1 , KFC 
FC( I )— 0 . 

FS( I )— 0. 

110  CONTINUE 
CALL  FOCFT 
DO  120  1-1  .  N 
II-  1+1 

DO  115  J-2.KFC 
J1«  J-1 

AW( J1 )-AA( I , J) 

BW(J1)«B3(I.J) 

CW( J 1 )*. 5» (  C(I,J)  +  C( 1 1 ,J)  ) 

DW(J1 )-. 5*(  0(1 ,J)  +  D( 1 1 ,J)  ) 

EW( J 1 )- . 5* (  A( I , J )  +  A( 1 1 , J )  ) 

FHK(J1)-.5*(  B(I.J)  +  B(  1 1  ,  J  )  ) 

115  CONTINUE 

AO-AA(I.1)«.5 

C0-.25»(  C(I,1)  +  C( 1 1 ,1)  ) 

AOO- . 25 • (  A ( I , 1 )  +  A ( 1 1 . 1 )  ) 

CALL  COED( AO ,C0, AOO ,U0) 

UCC( I , 1 )-2 . »U0 


&  1  OL. 
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DO  116  J-2.KFC 
J1-J-1 

UCC(I.J)«AA2(J1) 

UCS(  I .  J  )-8B2( J 1 ) 

VCC(I.J)*CC2(J1) 

VCS(I.J)-0D2(J1) 

116  CONTINUE 
120  CONTINUE 

DO  130  J-1 .N 

DO  129  1-2. IM2  , 

FC( I )»FC( I)  +  (VCS(J.n-UCC(J.l))  *CCP(I,J) 
FS(l)-FS(I)  +  (UCS(J.I)+VCC(J.O)  »CCP(I,J) 

129  CONTINUE 

FC(KFC)— FC(KFC)  -  UCC(J ,KFC)»CCP(KFC. J) 

FC(1 )— FC(1 )  +  UCC(J,1)»CCP(1 .J) 

130  CONTINUE 

DO  131  J-1.JDIM 
DO  131  I-1.IDIM 

04(1. J)-CCP(I ,J)/.4+.5*(U(I.J)**2+V(I.J)**2) 

131  CONTINUE 

DO  140  1-2. IM2 

FC(  I  )-FC(  I )  +VSC*B8(1,n 

FS( I )— FS(I )  +  VSC»AA( 1 , I ) 

140  CONTINUE 

DO  150  J-2.KFC 
FC(J)-FC(J)*SGMA(J) 

FS(J)-FS(J)*SGMA(J) 

150  CONTINUE 

DO  160  1-1 .IM 
CP(I)-1.-FC(1) 

DO  159  L-2.IM2  ,  ,  v  . 

159  CP(I)-CP(I)  +  2. *(  FC(L).CS(L.I)-FS(L).SN(L.I)  ) 
CP(I)-CP(I)  +  FC(KFC).CS(KFC.I) 

160  CONTINUE 
CP(I»iN-1)^:P(1) 

RETURN 

END 


*  »'■  V* 
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SUBROUTINE  COED (AO, CO, AOO.UO) 


MULTIPLICATION  OF  TWO  FOURIER  SERIES 
CALLS  00  NONE 
CALLED  BY  00  CPVAL 


PARAMETER  (IDIM-80) 

PARAMETER  (I CP-4-0) 

C0M40N/GRD/IM, IM2.KFC, JR,N 

C0»O40N/CPW/AW(  IDIM)  ,BW( IDIM) ,CW( ICP) ,DW( ICP) ,  EW(  ICP)  ,  FW(ICP) , 
1  AA2( ICP) ,BB2( ICP) ,CC2( ICP) , DD2( ICP) 

NM1—  IM2-1 
UO-  CO*AO 
00  100  I-1.IM2 

100  UO-  UO  +  . 5*(AW( I ) *CW( I )+BW( I ) »0W( I ) ) 

DO  110  1-1 , IM2 
II-  I+I 

AA2( I )-  AW( II) *CW(  I )  +  BW( II) «DW( I ) 

BB2( I )—  BW( II) *CW( I )  -  AW( 1 1 ) *DW( I ) 

CC2( I)-  AW(II)*EW(I)  +  8W(II)»FW(I) 

DD2( I )—  BW(II)*EW(I)  -  AW( 1 1 )  *FW( I ) 

110  CONTINUE 

DO  120  J-1 , NM1 
K4—  IM2-J 
DO  119  K— 1 ,J 
K1-  IM2+1-K 
K2—  K1— K4 
K3—  K1+K4 
W1-  AW(K2)+AW(K3) 

W2—  8WK2)+BW(K3) 

VK3-  AW(K2)-AW(K3) 

W4-  BW(K3)-BW(K2) 

AA2(K4)—  AA2(K4)  +  W1*CW(K1)  +  W2*DW(K1) 

BB2(K4)—  BB2(K4)  +  W3«DW(K1)  +  W4*CW(K1) 

CC2(K4)«  CC2IK4)  +  W1*EW(K1)  +  W2*FW(K1) 

DD2(K4)—  DD2(K4)  +  W3*FW(K1 )  +  W4*EW(K1 ) 

119  CONTINUE 

120  CONTINUE 

DO  130  J-1 , NM1 
K4—  J+1 
DO  129  K— 1 . J 
K2-  K4— K 
K3—  K4+K 

W1—  AW(K2)+AW(K3) 

W2—  BW(K3)-BW(K2) 

W3—  AW(K2)-AW(K3) 

W4—  BW(K2)+BW(K3) 

AA2(K4)»  AA2(K4)  +  W1*CW(K)  +  W2*0W(K) 

B82(K4)«  BB2(K4)  +  W3«DW(K)  +  W4*CW(K) 

CC2(K4)»  CC2(K4)  +  W1*EW(K)  +  W2*FW(K) 

DD2(K4)-  DD2(K4)  +  W3*FW(K)  +  W4*EW(K) 

129  CONTINUE 

130  CONTINUE 

00  140  1-1 ,IM2 

AA2( I )—  . 5*AA2( I )  +  AO»CW(I)  +  CO*AW( I ) 

BB2(  I )-  .  5*BB2(  I )  +  AO*OW(I)  +  CO*BV»(I) 

CC2( I )»  . 5*CC2( I )  +  AO*EW( I )  +  AOO*AW(I) 

DD2( I )-  .5*002(1)  +  AO*FW(I)  +  AOO*BW(I) 


U0  CONTINUE 

AA2(IM2)«  2.*AA2(IM2) 
CC2(IM2)-  2. «CC2( IM2) 
RETURN 


o  o  o  o  o  ooooooooooooooooo 


PROGRAM  LOADS 


LOADS  COMPUTATIONS  -  DISSPLA  VERSION 
LAST  REVISION  4-7-87 

PRINCIPAL  INVESTIGATOR  :  DR.  J.C.  WU 
AUTHORS  :  MIKE  PATTERSON,  ISHMAEL  TUNCER 
GEORGIA  INSTITUTE  OF  TECHNOLOGY 
(404)  894-3028 

TAPE3  :  INPUT  FROM  GEOM 

TAPE4  :  INPUT  FROM  ZONST 

TAPE5  :  GENERAL  INPUT 

TAPE6  :  GENERAL  OUTPUT 

TAPE10  :  OUTPUT  TO  PLOT3 


PARAMETER  ( IDU4-80 ,  JDIkA-60) 

PARAMETER  ( IP1-81 , JP1-61 ) 

DIMENSION  WB(IP1),CP(IP1),G(IP1) 

DIMENSION  XB( IDIM) , YB( IDIM) , XT( IDIM) , YT( IDIM) . FUN( IP1 ) 
COMKON/L1/DTET , IM, IM1 

REWIND  3 
REWIND  4 
REWIND  5 
REWIND  6 
REWIND  10 

READ(5 . »)  NTOUT.NTS.NTL 
RFAD(3)  AL.RE.OTET.XB, YB.XT, YT 

IM-IDIM 

JR-JDIM 

IM1-IP1 

JR1-JP1 

PI-4.«ATAN(1 .) 

C.  READ  SURFACE  VORTICITY  4  PRESSURE  COEFFICIENT 

100  READ( 4. END-99)  NT .NTPL.T ,RF,WB ,CP ,OMG ,OMGO , ALP 
IF(NT.LT.NTS)  GO  TO  100 
IF(NT.GT.NTL)  GO  TO  99 
ALPD-ALP* 1 80 . /PI 
WRITE(6 , 60)  NT.T.ALPD 
COS A-COS (ALP) 

SINA— SIN(ALP) 

C.  COMPUTE  TOTAL  CP  AND  OUTPUT  VALUES 


DO  110  1-1 , IM 

110  FUN( I )-OMGD* (XB( I ) «YT ( I )-YB( I ) - XT ( I ) ) 

FUN(IM1 )-FUN(1 ) 

CALL  INTER(FUN,G, 1 . IM1 ) 

DO  120  1-2, IM1 

WW-0MG»0MG»(XB(I)»XB(I)+YB(I)*YB(I)-XB(1)»XB(1 )) 

120  CP( I )-CP( I )-G( I )+WW 

I F(MOO(NT , NTOUT) . EQ. 0)  WRITE(6,61)  (CP( I ) , 1=1 . IM) 

C.. COMPUTE  TANGENTIAL  FORCE  COEFFICIENT  FROM  VISCOUS  EFFECT 


'Sgi'a'aWaBWro 
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DO  130  1-1 . IM1 
130  W8( I )-W8( I )— 2 . »OMG 
DO  140  1  =  1  , 1M 
140  FUN(I)-WB(I)»XT(I) 

FUN(1M1)«FUN(1) 

CALL  INTER( FUN ,G,  1 ,  IM1  ) 

CTV-2 . 0»G( 1M1 )/RE 

C.. COMPUTE  NORMAL  FORCE  COEFFICIENT  FROM  VISCOUS  EFFECT 

DO  150  1-1 . IM 
150  FUN(I)«W8(I)»YT(I) 

FUN(IMI)— FUN(1) 

CALL  INTER(FUN,G, 1 ,  IM1 ) 

CNV— 2 . 0*G( IM1 )/RE 

C.. COMPUTE  MOMENT  COEFFICIENT  FROM  VISCOUS  EFFECT 
DO  160  1-1  ,IM 

160  FUN( I )-WB( I ) • (XB(I )«YT ( I )— Y8( I ) -XT ( I ) ) 

FUN( IM1 )— FUN( 1 ) 

CALL  INTER(FUN.G.I.IMI) 

CMV— 2 .0*G( 1M1 )/(RE»AL) 

WRI TE(6, 66)  CNV.CTV.CMV 

C.. COMPUTE  TANGENTIAL  FORCE  COEFFICIENT  FROM  PRESSURE  DISTRIBUTION 

DO  170  1-1 .IM 
170  FUN( I )-CP( I ) »YT ( I ) 

FUN(IM1)-CP(IM1)«YT(1) 

CALL  INTER(FUN.G.I.IMI) 

CTP— G(  IM1  )/AL 

C. . COMPUTE  NORMAL  FORCE  COEFFICIENT  FROM  PRESSURE  DISTRIBUTION 


DO  180  1-1 .IM 
180  FUN(I)-CP(I)»XT(I) 

FUN( IM1 )— CP( IM1 )»XT(1  ) 

CALL  INTER(FUN.G.I.IMI) 

CNP-G(IM1)/AL 

C.. COMPUTE  MOMENT  COEFFICIENT  FROM  PRESSURE  DISTRIBUTION 
DO  190  1-1 .IM 

190  FUN(I)-CP(I)*(XB(I)*XT(I)+YB(I)*YT(I)) 

FUN( IM1 )-CP( IM1 )«(XB(1 )*XT( 1 )+YB( 1 ) »YT( 1 ) ) 

CALL  INTER(FUN.G.I.IMI) 

CMP— G( IM1)/(AL«AL) 

C .  WRITE  PRESSURE  COEFFICIENTS  AND  COMPUTE  TOTAL  NORMAL,  TANGENTIAL 
C  AND  MOMENT  COEFFICIENTS 

WR I TE(6 . 64)  CNP.CTP.CMP 
CN-CNV+CNP 
CT-CTV+CTP 
CM— CMV+CMP 


C.. COMPUTE  LIFT  AND  DRAG  COEFFICIENTS 

CLP-CNP«COSA-CTP*SINA 
CDP-C  NP • S I N A+C T P  «COS A 
C L=CN*COSA-CT • S I NA 
CD-CN* S I NA+CT «COSA 


ooooooooo 


WRITE(6 . 67)  CLP, COP, CMP 
WRITE(6 , 69)  CL, CD, CM 

C.. WRITE  DATA  FOR  PLOT 3  PROGRAM  AND  RETURN  TO  BEGINNING  OF  PROGRAM 

WRITE(10)  NT ,NTPL , T , ALPD , RF , RE , CLP ,CDP ,CMP ,CP 
GO  TO  100 

C . . LAST  OF  DATA  HAS  BEEN  READ 
99  WRITE(6 , 63) 

60  FORMAT ( 1 H  ,40(1W-),/-  NT  -' , 14, ’ .  T  -• ,F7.4 

+  AA  -•  .F8.4/1X, 40(11+-)) 

61  FORMAT (// . ‘  CP  VALUES.  1-1 . IM  : '//( 1 0E13 . 6) ) 

63  FORMAT (//’  - :  END  OF  FILE  : - '//) 

64  FORMAT ( ’  VALUES  OF  CNP.CTP.CMP  - ’.5X.3E16.8) 

66  FORMAT (//'  VALUES  OF  CNV.CTV.CMV  - - — ',5X,3E16  8) 

67  FORMAT ('  PRESS  LOADS  CLP ,C0P , CMP  — 1 , 5X , 3E1 6 , 8) 

69  FORMAT ( '  TOTAL  LOADS  CL, CD. CM  - ‘ ,5X,3E16,8) 

STOP 
END 

SUBROUTINE  INTER(A,AI . IS, IL) 


NUMERICAL  QUADRATURE  BY  HIGH  ORDER  FORMULAS 
CALLS  :  NONE 
CALLED  BY  :  LOADS 


DIMENSION  A( IM1 ) . AI ( 1M1 ) 

C0M40N/L1/DTET .  IM,  IM1 

INC-1 

CI-DTET/24. 

IF(IL.LT.IS)  THEN 
INC— INC 
Cl— Cl 
END  IF 
AI ( IS)— 0 . 

DO  100  1-IS+INC. IL, INC 
IM2-I-2* INC 
IP1-I+INC 

I F( IM2 . LT . 1 .OR. IM2.GT.IM1)  IM2-ABS( IM1-1-IM2) 

IF(IP1 . LT. 1 .OR. IP1 .GT. IM1 )  IP1-ABS( IM1-1-IP1 ) 
AI(I)«AI(I-INC)+CI.(-A(IM2)+13.»(A(I-INC)+A(I))-A(IP1 )) 
100  CONTINUE 
RETURN 
END 


PROGRAM  PLOT  1 


GEOMETRY  PLOTTING  PROGRAM 
LAST  REVISION  4-1-87 

PRINCIPAL  INVESTIGATOR  :  OR.  J.C.  WU 
AUTHORS  :  MIKE  PATTERSON,  ISHMAEL  TUNCER 
GEORGIA  INSTITUTE  OF  TECHNOLOGY 
(404)  894-3028 


TAPE1 

TAPE5 

TAPE11 


INPUT  FROM  GEOMETRY  PROGRAM 
GENERAL  INPUT 
OUTPUT  FOR  PLOTTER 

AXES.  GRID.  WAKE.  ZONES 


PARAMETER  (IDIM-80 ,  JOIIA-60) 

PARAMETER  ( IP1-81 , JP1-61 ) 

PARAMETER  (I NOR 1-20) 

C0M40N/GRD/X ( I D I M , J P 1 ) . Y ( I D I M , J P 1 ) 

C0M40N/ZN/AL.JVL, J8L.JFL.IV1 .IB1 . IB2 , I V2 . I LINE(6) , SCALE 
C0M40N/WK/IN0R(IN0R1 ,JDIM) . IWK( JDIM. JDIM) 

C0M40N/P  AR  AM/ 1 M ,  J  R 
DATA  IRES/2/ 

REWIND(I) 

REWIND(5) 

READ ( 1 )  AL.X.Y, I NOR , IWK 
READ(5, *)  IV1 . IB1 . IB2, IV2 
REA0(5..)  IPLOT1 . IPLOT2. IPL0T3, IPL0T4 
READ(5.«)  JVL.JBL.JFL, SCALE. IDEV 

IM-IDIM 

JR-ODIM 

I LINE( 1 )— 1 

ILINE(2)-IV1 

ILINE(3)-IB1 

ILINE(4)— IB2 

ILINE(5)-IV2 

ILINE(6)-IM 

SET  DEVICE 

CALL  D IP( 1 5 . 'PLOT  1 .DIP' .9) 

1F(IDEV.EQ.0)  THEN 
CALL  EPIC(300. .0,8  25,10.75.11) 

ELSE  I F( IDEV  EQ  1  )  THEN 
CALL  PCLCfcP 

CALL  CALCInf,(IBUF,512.11) 

ELSE  I F( [DEV  EQ.2)  THEN 
CALL  PVRSTC 

CALL  VRSTEC( IBUF ,512,11) 

ELSE  I F( IDEV  EQ  3)  THEN 
CALL  P4115 
CALL  TK4 1 15( IRES) 

END  IF 


C  CONSTRUCT  EACH  PLOT  SPECIFIED 


o  o  o  o  o  o 


c 

AIRFOIL  AND 

FLOW  ZONES 

V 

c 

CALLS 

AXES.  GRID 

c 

CALLED  BY  : 

PLOT1 

i 

c 

PARAMETER  1 

;iDIM-80, JOIM-60) 

-1 

PARAMETER  ( 

’  IP1-81 , JP1— 61 ) 

C0M40N/GRD/X(IDIM,  JP1  )  .Y(IOIM.  JP1  ) 

C0MM0N/ZN/AL.JVL.J8L.JFL.IV1 .181 ,182, IV2.ILINE(6). SCALE 
COMMON/P ARAM/ IM . JR 
DIMENSION  XX ( IP1),YY(IP1) 

CALL  AXES(JFL.AL/4. , ' FLOW  ZONES$’) 

CALL  THKFRM(.01) 

CALL  FRAME 
CALL  GRID(0, IM.5) 

CALL  HEIGHT (.1) 

C.  PLOT  AND  NUMBER  THE  DEMARCATION  LINES 

DO  110  1-1,6 
DO  109  J-1 , JFL 


xx(j)-x(iline(i),j: 

YY(J)-Y(ILINE(I)  ,  J) 


109  CONTINUE 

CALL  CURVE ! XX ,YY , JFL.0) 

CALL  RLINT(ILINE(I) ,XX(JFL) .  Y( I LINE( I ) , JFL+1  ) ) 

1 10  CONTINUE 

CALL  HEIGHT( . 14) 

CALL  ENOPL(0) 

RETURN 

END 

SUBROUTINE  WAKE 


AIRFOIL  AND  WAKE  GRID 
CALLS  :  AXES.  GRID 


riT«-inn^ 


c 

PARAMETER  ( IDIM-80 , JDIM-60) 

PARAMETER  ( I P 1 =8 1 . JP1-61 ) 

PARAMETER  (I NOR 1-20) 

COMMON/GRD/X( IDIM.JP1) .Y(IDIM,JP1 ) 

COMHON/ZN/AL.  JVL,  JBL,  JFL.  I VI  ,  IB1 . I B2 . IV2 , I L1NE(6) , SCALE 
COMi»N/WK/INOR(INOR1  .  JOIM)  .  IWK(JDIM,  JDIM) 

COMADN/P  ARAM/ 1 M ,  J  R 
DIMENSION  XX(IP1).YY(IP1) 

C. .PLOT  AXES  AND  AIRFOIL 

CALL  AXES ( JVL. AL. ‘WAKE  GRIDS') 

CALL  GRID(0, IM,5) 

C . . PLOT  VORTICITY  GRID  UNDERNEATH  WAKE  GRID 
CALL  DOT 

CALL  THKCRV( . 001 ) 

CALL  GRID( IM/4+1 , I M/4+1 , JR) 

CALL  RESET ('DOT') 

CALL  RESET ( 'THKCRV' ) 

C . . PLOT  STREAMWISE  LINES 

DO  1 10  Jt+-1  .  JR-1 
J  J— 0 

DO  109  J-JN+1 , JR 
JJ-JJ+1 

XX(  JJ )*X( IWK(JN.J).J) 

YY(JJ)-Y(IWK(JN,J) . J) 

109  CONTINUE 

CALL  CURVE(XX . YY , J J , 0) 

110  CONTINUE 

C  PLOT  NORMAL  LINES 

DO  120  1—2. IM/4 
DO  119  J-1 .JR 

XX( J )— X( INORf I . J) , J ) 

YY(J)-Y(INOR(I. J). J) 

119  CONTINUE 

CALL  CURVE(XX.YY. JR.0) 

120  CONTINUE 
CALL  ENDPL(0) 

RETURN 

END 

SUBROUTINE  AXES( JMAX .OFFSET . LABEL) 

C 


ooooooooo 


dll' 


COMMON/ZN/AL. JVL. JBL. JFL. IV1 . 181 . IB2, 1V2, I LINE( 6) , SCALE 
COMMON/PAR AM/ 1 M ,  J  R 
CHARACTER *20  LABEL 


C . . PAGE  LAYOUT 


PAGEX-8 . 25 
PAGEY-10.75 
XAXIS-7.25 
YAXIS-7.25 
CALL  BLOWUP(SCALE) 

CALL  PAGE(PAGEX.PAGEY) 

CALL  AREA2D ( XAX I S , YAX I S ) 

CALL  HEAD1N(LABEL, 100.1 .25,1) 


C.  COMPUTE  GRID  AXIS  LENGTHS  AND  SCALES  BASED  ON 
C  MAXIMUM  VALUES  TO  BE  PLOTTED 


GRIDMAX-0. 

DO  110  1-1 . IM 
DO  109  J~1 . JMAX 


I F  ( A8S ( X ( I , J ) )  ,GT.  GRIDMAX)  GRIDMAX»ABS(X( I . J ) ) 
IF(ABS(Y(I ,J))  GT.  GRIDMAX)  GRIDMAX-ABS(Y( I . J ) ) 


109  CONTINUE 

110  CONTINUE 


C.  LOCATE  THE  ORIGIN  OF  THE  PLOT  AND  DRAW  AXES 


XOR I G— GR I DMAX+OFFSET 
YOR I G— GRIDMAX 
XMAX-GR I  DMAX-fOFFSET 
YMAX-GRIDMAX 
XSTP-1 . 

YSTP— 1 . 

CALL  GRAF (XOR I G . XSTP . XMAX . YORIC . YSTP . YMAX) 

RETURN 

END 

SUBROUTINE  GRIO( IMAX1 . IMAX2, JMAX) 


GRID  IN  PHYSICAL  PLANE 


CALLED  BY  :  PLOT 1 .  ZONES.  WAKE 


PARAMETER  ( IDII*«80  .  JDlkA-60) 
PARAMETER  ( IP1-81 , JP1-61 ) 


COMWON/GRD/X( IOIM, JP1 ) . Y( ID IM. JP1 ) 
COMMON/P ARAM/ 1 M . J  R 
DIMENSION  XX( IP1 ) . YY( IP1 ) 


C.  PLOT  THE  RADIAL  LINES 


DO  110  1-1 . I MAXI 
DO  109  J-1 .JMAX 
XX(J)-X(I ,J) 

YY( J )— Y( I , J ) 

109  CONTINUE 

CALL  CURVE(XX.YY. JMAX.0) 

1 10  CONTINUE 


177 
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C •  PLOT  AZIMUTHAL  LINES 


II-IMAX2 
DO  120  J— 1 , JMAX 
DO  119  1-1 . IMAX2 
XX(I)-X(I.J) 
YY(I)-Y(I.J) 

119  CONTINUE 

IF( IMAX2. EQ. 1M)  THEN 
XX( IMAX2+1 )«XX( 1 ) 
YY(IMAX2+1)-YY(1) 

1 1— IMAX2+1 
END  IF 

CALL  CURVE ( XX , YY , 1 1 , 0 ) 

120  CONTINUE 
RETURN 
ENO 
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PROGRAM  PLOT 2 


THIS  PROGRAM  ACCEPTS  20-ARRAYS  TO  MAKE  CONTOURS 
LAST  REVISION  5-25-86 


PRINCIPAL 

INVESTIGATOR  :  DR.  J.C.  WU 

AUTHOR 

MIKE  PATTERSON 

GEORGIA  INSTITUTE  OF  TECHNOLOGY 

(404)  894-3028 

TAPE5 

GENERAL  INPUT 

TAPE6 

GENERAL  OUTPUT 

TAPE9 

INPUT  FROM  ZONST  FOR  STREAMLINES  AND  CONTOURS 

TAPE  12 

OUTPUT  FOR  PLOTTER 

TAPE14  : 

INPUT  FROM  GEOM 

CALLS 

CONTG 

PARAMETER  ( ID  11^80 , JOIM-60) 

PARAMETER  ( IP1-81 , JP1-61 ) 

REAL  XX(IP1 ,JP1),YY(IP1,JP1),S1(IP1 . JP1 ) ,S2( IP1 . JP1 ) 
REAL  XB(IDIM) ,Y8(I0IM) ,CONT(100) 

REAL  XB1 ( IP1 ) ,  YB1 ( IP1 ) 

INTEGER  IBUF(51 2) 

COMMON/SUM/CSQ . GAMA .SIGMA, COSA , S I NA 
DATA  IRES/2/ 

DATA  PI/3.14159/ 

REWIND  5 
REWIND  6 
REWIND  9 
REWIND  14 

CALL  DIP(  15.  ’PLOT2.DIP' .  9  ) 

.  GENERAL  INPUT 

READ(5 , • )  VALMIN, VALMAX.NCON 

READ(5 . * )  VORMIN , VORMAX , NVOR 

READ(5 , • )  AH I . ICHR .NDIG .NSKIP . LINTP .HEAVY 

READ(5. •)  XI .X2.DX.XLEN 

READ(5 , • )  Y1 .Y2.DY.YLEN 

READ(5, •)  IDEV, IOPT . IPRINT, IPAR 

READ(5 , « )  JMAX. SCALE. NTS, NTL 

IF(NCON.GT.0)  READ(5 . • )  (CONT( I ) , 1-1 . NCON) 

IVMDIM 

IM1-IP1 

.  . INPUT  FROM  GEOM 

RCAD( 14)  CSO , GAMA , S I GMA , AL 
READ( 14)  XB.YB.XX.YY 

.  .SET  OEVICE 

I F( IDEV . EO . 0)  THEN 
CALL  EPIC(300. ,0.8.25. 10.75. 12) 

ELSE  IF(IDEV.EQ.I)  THEN 
CALL  PCLCMP 

CALL  CALCMP( IBUF. 51 2 .12) 


o  o  o  o  o  o  o 


ELSE  IF(ID£V.E0.2)  THEN 
CALL  PVRSTC 

CALL  VRSTEC( I8UF ,512,12) 

ELSE  I F ( IDEV . EQ . 3)  THEN 
CALL  P4115 
CALL  TK41 15( IRES) 

END  IF 

CALL  SETDEV(6,6) 

C. .INPUT  FROM  ZONST  —  OUTER  LOOP  —  NEW  PAGE  OF  PLOTS 
1000  CONTINUE 

READ (9. END-3000)  NT,NPL,T,ALP0,RF,RE,S2,S1 
IF(NT.LT.NTS)  GO  TO  1000 
IF(NT.GT.NTL)  GO  TO  3000 
WRITE(6, 10)  NT,NPL,T,ALPD,RF,RE 
IF(JMAX.GT.NPL)  THEN 

WRITE(6.  »)  '  JMAX  >  NPL  - EXPECT  GARBAGE - •’ 

STOP 
END  IF 

ALP— ALPD»PI/180 . 

SINA— SIN(ALP) 

COSA-COS(ALP) 

I SUB-0 
I NAME- I OPT 
CALL  NOBRDR 
CALL  NOCHEK 
CALL  GRACE(0. ) 

C.  INNER  LOOP  —  SUBPLOTS 

2000  CONTINUE 

ISUB-ISUB+1 

IF(ISUB  EQ. 1 )  CALL  PHYSOR(1 .5.5.50) 

IF( ISUB. EQ. 2)  CALL  OREL(0. .-4.25) 

CALL  COMP LX 

CALL  BASALF( 'STAND') 

CALL  MIXALF( ’ L/CSTD’ ) 

CALL  MX3ALF( ' L/CGR ’ . 1H/) 

CALL  BLOWUP(SCALE) 

CALL  AREA2D(XLEN.YLEN) 

CALL  YAXANG (0.0) 

CALL  XTICKS(0) 

CALL  YTICKS(0) 

CALL  XNONUM 
CALL  YNONUM 

CALL  GRAF(X1 .DX.X2.Y1 .DY.Y2) 

CALL  THKFRM( . 02) 

CALL  FRAME 

C.  P LOT  PARAMETERS 

I F ( I PAR . EQ . 1 )  THEN 
WB  -  1 . 55 
HB  -  1 . 05 
X8LK  -  .05 
YBLK  -  .05 

CALL  BLTREC(XBLK , YBLK ,WB ,HB,0 . 0 , 0 . 01  5) 

CALL  BLKEY(ID) 

CALL  BLOFF(ID) 

XM  -  XBLK  +  0.10 


YM  -  YBLK  +0.05 
CALL  HEIGHT(0. 15) 

CALL  MESSAG('R(E  -)$  \  100 , XM, YM) 

XR-XM+- .  5 

CALL  R£ALN0(R£,-1 ,XR. YM) 

XM  -  XBLK  +  0.1 

YM  -  YBLK  +0.30 
CALL  H£IGHT(0. 15) 

CALL  MESSAG('(RF  -)$\  100.XM.YM) 

XR  -  XM  +  0.6 
CALL  R£ALN0(RF,3.XR.YM) 

XM  -  XBLK  +  0.1 

YM  -  YBLK  +0.55 
CALL  HEIGHT (0.15) 

CALL  MESSAG(’(T  -)$ * . 100 . XM. YM) 

XR  -  XM  +  0.6 
CALL  REALN0(T.3.XR.YM) 

XM  -  XBLK  +  0.1 

YM  -  YBLK  +  0.80 
CALL  HEIGHT (0.15) 

CALL  MESSAGC/A)  «$' .100.XM.YM) 

XR  -  XM  +  0.6 
CALL  REALNO( ALPO , 3 . XR , YM) 

CALL  BLON(ID) 

END  IF 

C. .DRAW  BLADE  AT  ANGLE  OF  ATTACK 
DO  140  1-1. IM 

X81 ( I )-XB( I ) «COSA+YB( I)*SINA 
YB 1 ( I )-Y8 ( I ) .COSA-XB ( I ) • S I NA 
140  CONTINUE 

XB1 ( IM1 )-XB1 ( 1 ) 

YB1  (  IM1 )— YB1  (  1 ) 

CALL  THKCRV(0 . 010) 

CALL  CURVEfXBI ,YB1 . IM1  .0) 

CALL  RESET ( ’ THKCRV ’ ) 

C.. DISPLAY  DATA.  CALL  CONTOURING  ROUTINE 

IF( IPRINT . EQ. 1 ) 

1  WRITE(6,20)((I,J ,XX(I ,J),YY(I,J),S1(I,J),S2(I,J), 

2  J— 1 . JMAX.2) , 1—1 . IM,4) 

IF(ISUB.EQ.I)  THEN 

CALL  CONTG(S1 .XX.YY.1 . IM1 .1 , JMAX , VALMIN , VALMAX , AH I .CONT.NCON. 
1  ICHR . NDIG.NSKIP , LINTP , HEAVY) 

ELSE 

CALL  CONTG(S2 , XX . YY , 1 . IM1 .1 . JMAX , VORMIN , VORMAX . AHI .CONT.NVOR. 
1  ICHR. NDIG.NSKIP, LINTP. HEAVY) 

END  IF 

C.. WRITE  APPROPRIATE  ] ITLE 

XM-XLEN-4.0 
YM-YLEN- . 25 
I F ( I NAME . EQ . -2 )  THEN 

CALL  MESSAGO  S ( TREAML I NES ) $  M00.  XM,  YM) 

ENOIF 

IF(INAME.EO.-I)  THEN 

CALL  MESSAG( ’V(ORTICITY)  C(ONTOURS) V . 1 00 . XM, YM) 
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ENDIF 

I F ( I NAME .  EQ . 0 )  THEN 

CALL  MESSAG ( ' 0 ( ENS I TY )  C(ONTOURS)$ 1 . 1 00 , XM. YM) 

ENDIF 

I F ( I  NAME . EO . 1 )  THEN 

CALL  MESSAG ( 'M(ACH)  N(UMBER)  C(ONTOURS)$' , 1 00 , XM, YM) 

ENDIF 

I F( I NAME . EO . 2 )  THEN 

CALL  MESSAG ( ’ V ( ELOC I TY )  M(AGNITUDE)  C(ONTOURS)$‘ . 100.XM, YM) 
ENDIF 

IF(INAME.EQ.3)  THEN 

CALL  MESSAG( ’S(KIN)  F(RICTION)  C(ONTOURS)$‘ , 100.XM.YM) 

ENDIF 


C. .FINISH  SUBPLOTS  OR  START  NEW  PAGE 


CALL  ENOGR(ISUB) 

INAME-INAME+1 
I F ( I SUB . EQ . 1 )  GO  TO  2000 
CALL  ENOPL(0) 

GO  TO  1000 

C.. PLOTTING  FINISHED 

3000  CALL  DONEPL 
STOP 

10  FORMAT! '  1 . ’NT.NPL.T.ALPD.RF.RE' .2I4.4F13.3) 

20  FORMAT(2(5X,2I5 , 4F10 . 5) ) 

END 

SUBROUTINE  CONTG(A,XA, YA, IL. 1U. JL, JU.VALMIN.VALMAX. 

1  AHI.CONT.NCON. ICHR . NDIG .NSKIP . LINTP .HEAVY) 

C 

C  . . . . . . * . . . 

C  THIS  SUBROUTINE  PLOTS  CONTOURS  OF  MATRIX  "A”  FOR  GENERAL  (X.Y) 
C  POSITIONS  OF  THE  ELEMENTS 

C 

C  CALLS  :  XFUN 

C 

C  CALLED  BY  :  PLOT2 

C  . . . . . . . . 

c 

PARAMETER  ( IP1-81 , JP1-61 ) 

DIMENSION  A(IP1,JP1),XA(IP1 , JP1 ) . YA( IP1 , JP1 ) ,CONT( 100) 
DIMENSION  B(4) . X(4) ,  Y(4) 

COfcMON/COORD/XX(3) . YY(3) ,XX2(3) ,YY2(3) 

REAL  HT(4) 

LOGICAL  HVY 

DATA  HT/.08. .10. . 12. .14/ 

DATA  PI/3.14159/ 

IF( ICHR . GT . 0)  CALL  HEIGHT(HT( ICHR) ) 

I2-IU-1 

J2-JU-1 

NSKP-NSKIP 

I F( LINTP . EO  3)  CALL  CHNDOT 
I F( LINTP . EQ . 4)  CALL  CHNDSH 
I F( LINTP . EQ . 2)  CALL  DASH 
I F( LINTP. EQ. 1 )  CALL  DOT 
IF (NSKIP. LE . 0)  NSKP-10 
HVY-. FALSE. 

IF(VALMIN . LT . VALMAX )  THEN 
VALMX-VALMAX 
VALMN-VALMIN 

ELSE 
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VALAW-1 £30 
VALMX— 1E30 
00  170  I-IL.IU 
DO  170  J-JL.JU 
IF  (A(I.J).LT.AHI)  THEN 

VALMX-AMAX1 (A( I . J ) . VALMX) 

VALMN-AMIN1 (A(I ,J) , VALMN) 

END  IF 

170  CONTINUE 

END  IF 

kWKiT'AKI 

I F(NC0N. LT . 0)  THEN 
N  -  -NOON 
CONMIN  -  VALMIN 

CON I NO  -  (VALMAX  -  VALMIN)/FL0AT(N-1 ) 

ENDIF 

C. .DO  LOOP  FOR  CONTOUR  LINES 

DO  500  NCN-1 , N 
I F ( NOON . LT . 0 )  THEN 

CONTUR-CONM I N+ ( NON- 1 )»CONINC 
ELSE 

CONTUR-CONT(NCN) 

ENDIF 

I F(HEAVY . NE . 0 . ) 

1  HVY«A8S(AINT(ABS(CONTUR/HEAVY)+0.5)-A8S(CONTUR/HEAVY) ) . LT  0  001 
IF(HVY)  CALL  THKCRV(.01) 

NSK-NSKP 

C. SEARCH  MATRIX  "A"  FOR  CONTOURS 

DO  490  I-IL.I2 
DO  490  J-JL.J2 

B8-AMAX1(A(I.J).A(I ,J+1),A(I+1 ,J),A(I+1 ,J+1)) 

IF  (BB-AHI )  290.490.490 
290  IF  (CONTUR-BB)  300.490.490 

300  IF  (CONTUR— AMIN1 (A( I ,J) ,A(I ,J+1 ) ,A(I+1 ,J) , A( 1+1 , J+1 ) ) )  490.310.310 
310  B(1)».25*(A(I ,J)+A(I+1 . J )+A( I , J+1 )+A( 1+1 ,J+1)) 

B(4)«B( 1 ) 

X( 1 )- . 25«(XA( I . J )+XA( 1+1 . J )+XA( I . J+1 )+XA( 1+1 . J+1 ) ) 

X(4)«X(1) 

Y(1)-.25*(YA(I ,J)+YA(I+1 . J )+YA( I , J+1 )+YA( 1+1 . J+1 ) ) 

Y(4)-Y(1) 

C.. SEARCH  MATRIX  SUB-CELL  TRIANGLES  FOR  CONTOURS 

DO  480  K-1  .4 
NP-0 

IF  (K.LE.1)  THEN 
B(2)-A(I+1 .J) 

X(2)-XA(I+1 , J) 

Y(2)-YA(I+1 . J) 

ELSE 

8(2)-B(3) 

X(2)-X(3 

Y(2)-Y(3) 

ENDIF 

IlnM+K/3 

JlnHJ+1-I  ABS(5-2*K)/3 
B(3)-A( IM. JM) 

X(3)-XA(IM. JM) 

Y(3)-YA(IM. JM) 


C.. DETERMINE  CONTOUR  INTERSECTIONS 


DO  430  M-1  ,3 

IF  (CONTUR— AMIN1  (B(M)  . 8(M+1 ) )  )  4-30,380.380 
380  IF  (CONTUR-AMAX 1 (8(M),B(M+1)))  390,390,430 
390  NP-NP+1 

BB— 8(M+1 )— B(M) 

IF  (ABS(BB) . LE. 1 . E-15)  THEN 
D-0.5 

ELSE 

D-(C0NTUR-8(M))/8B 
END  IF 

XX(NP)-X(M)+0*(X(Mf1)-X(M)) 

YY(NP)-Y(M)+0»(Y(Mf1)-Y(M)) 

430  CONTINUE 

IF(NP.GT.I)  THEN 
I F( ICHR. LT.0)  THEN 

I F (CONTUR . GT . 0 . )  CALL  RESET( 'DASH ’ ) 

I F (CONTUR .  LE . 0 . )  CALL  DASH 
END  IF 

C. .TRANSFORM  CONTOUR  TO  PHYSICAL  PLANE  AND  DRAW 

CALL  XFUN(NP) 

CALL  CURVE(XX2 , YY2 ,NP , 0) 

C. . LABEL  CONTOURS 

I F ( ICHR .GT . 0)  THEN 
NSK-NSK+1 

IF(NSK.GE.NSKP)  THEN 
NSK— 1 

XS-.5*(XX2(1)+XX2(2)) 

DXS— XX2 ( 2 )— XX2 ( t ) 

YS-.5*(YY2(1)+YY2(2)) 

DYS-YY2(2)-YY2(1) 

ANG- 1 80 .  *  AT  AN2  (  DXS ,  DYS  )  /P I 
CALL  ANGLE(ANG) 

CALL  REA LNO( CONTUR , NO I G , XS , YS ) 

END  IF 
END  IF 
END  IF 

480  CONTINUE 
490  CONTINUE 

IF(HEAVY.NE. 0)  CALL  RESET ( 1 THKCRV ’ ) 

500  CONTINUE 
RETURN 
END 

SUBROUTINE  XFUN(NP) 

C 

C  . . 

C  COMPUTES  CARTISIAN  COORDINATES  (XX, YY)  IN  PHYSICAL  PLANE 
C  OF  A  SPECIFIC  POINT  (X.Y)  IN  COMPUTATIONAL  PLANE 
C 

C  CALLS  :  NONE 

C 

C  CALLED  BY  :  PLOT2 ,  CONTG 

C  . . . . . . MM 

C 

COMMON/SUM/CSQ , GAMA , S I GMA , COSA . S I N A 
COfc*40N/COORD/X  (  3  )  ,Y(3)  ,XX(3)  .  YY(3) 

COMPLEX  Z1 ,Z 

DO  100  1-1 , NP 
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Z1-GAMA+CMPLX(X( I ) , Y( I ) ) 

Z-Z1+CSQ/Z1+SIGMA 

XI-REAL(Z) 

Y 1  »A  l MAG ( Z ) 

XX ( I )»X1 .C0SA+Y1 *SINA 
YY ( I ) -Y 1 .COSA-X 1  • S I NA 
100  CONTINUE 
RETURN 
END 


PROGRAM  PLOT 3 


c 

c 

LOADS  PLOTTING  PROGRAM 

c 

LAST  REVISION  6-25-87 

c 

c 

PRINCIPAL 

INVESTIGATOR  :  DR.  J.C.  WU 

c 

AUTHOR 

MIKE  PATTERSON 

c 

GEORGIA  INSTITUTE  OF  TECHNOLOGY 

c 

(404)  894-3028 

c 

c 

TAPES 

GENERAL  INPUT 

c 

TAPES 

GENERAL  OUTPUT 

c 

TAPE10 

INPUT  FROM  LOADS  PROGRAM 

c 

TAPE 13 

OUTPUT  TO  PLOTTER 

c 

TAPE14 

INPUT  FROM  GEOM 

c 

c 

c 

CALLS 

NONE 

PARAMETER  (ID1M-80, JDIM-60) 

PARAMETER  (IP1-81 , JP1-61) 

PARAMETER  ( IU-40. I L-42 , 1UAX-2000) 

REAL  XB( IDIM) , XU( IU) ,XL(IL) 

REAL  CP( IP1 ) ,CPU( IU) ,CPL( I L) 

REAL  T ( I MAX ) .CL(IMAX) .CD(IMAX) .CM(IMAX) , ANG ( I  MAX ) 

REAL  XEXP(IU).  CPEXP(IU.JDIM) .  CLEXP(IU).  CDEXP(IU). 

>  CMEXP (  I L' '  ,  ALPHA(IU) 

INTEGER  IBUF(512) ,NT( I MAX) 

DATA  IRES/2/ 

C  REWIND  5 

C  REWIND  6 

C  REWIND  10 

C  REWIND  14 

C . . READ  INPUT  PARAMETERS 

READ(5 . • )  SCALE. SCALE2 
READ(5, •)  XI .X2.DX.XLEN 
READ(5 . • )  Y1 .Y2.DY.YLEN 
READ(5,  •)  I  OPT ,  I  DEV .  I  PAR ,  NTS ,  NT  I, 

C. . INPUT  FROM  GEOM 

READ(14)  DUMM1 . DUMM2 , 0UfcM3 . AL 
READ( 14)  XB 

C.. INPUT  FROM  EXPERIMENTAL  RESULTS 

C  READ(24 , • )  NPARTS 

C  READ( 24 , 900)  (  XEXP( I ) , 1=1 . 1 6  ) 

C  DO  J-1, NPARTS 

C  READ(24 , 900)  ALPHA ( J ) ,  CLEXP(J).  CDEXP(J),  CMEXP(J) 

C  READ(24 . 900)  (  CPEXP( I . J ) , 1-1 .  1 6  ) 

C  DO  1-1 . 16 

C  CPEXP(I.J)  -  -CPEXP(I.J) 

C  ENODO 

C  ENODO 

C  I  EXP  -  0 

900  FORMAT (  5(  1X.E14  7  )  ) 
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C. .SET  DEVICE 

CALL  DIP(  15.  'PL0T3.DIP1 .  9  ) 


IF(IDEV.EO.O)  THEN 
CALL  EPIC(300. ,0.8.25,10.75,13) 
ELSE  IF(IDEV.EQ.I)  THEN 
CALL  PCLCMP 

CALL  CALCMP( IBUF ,512.13) 

ELSE  IF{ IDEV .£0.2)  THEN 
CALL  PVRSTC 

CALL  VRST£C( I8UF.512 , 13) 

ELSE  IF(IDEV.E0.3)  THEN 
CALL  P4115 
CALL  TK41 15( IRES) 

ENDIF 

CALL  SETDEV(6,6) 

IV^IDIM 


. . CP  PLOT  LOOP. 

I F( IOPT . EQ. 0)  THEN 
INPUT  FROM  LOADS 


1000  CONTINUE 

READ(10. END-2000)  NT1 .NTPL.T1 , ALPD.RF.RE.CL1 ,CD1 ,CM1 .CP 
IF(NT1 . LT.NTS)  GO  TO  1000 
IF(NT1  GT.NTL)  GO  TO  2000 
IF(MOO(NT1 .NTPL) .NE.0)  GO  TO  1000 
WRITE(6.10)  NT1.T1.ALP0 
10  FORMAT  ( '  •,'NT.T.ALPD  :M5.2F8.3) 


I EXP  -  I EXP  +  1 

..SET  UP  ALPHA-NUMERICS  ANO  AXES  NAMES 


CALL  NOBRDR 
CALL  COMP LX 
CALL  BASALFC STAND1) 
CALL  MIXALF( ' L/CSTD 1 ) 
CALL  MX3ALFI 1 L/CGR 1 . 1 H\) 
CALL  BLOWUP(SCALE) 


.  .  LABEL  AXES 


CALL  XNAMEf ’X/C$' .100) 
CALL  YNAME('-C(P)$\100) 


.DIVIDE  DATA  INTO  UPPER  AND  LOWER  ARRAYS 


DO  100  1-1 . IU 

XU( I )— (XB( I )+AL/4 . )/AL 
CPU(I)— CP(I) 

100  CONTINUE 

IF (  I  TEST . EO. 0  )THEN 
WRITE(24, • )  (  XU(I).I-IU.1,-1  ) 
ITEST-1 
ENDIF 

WRITE(24, • )  ALPD 

WRITE(24. •)  (  CPU(I).I»IU,1.-1  ) 


I M2- I M/2-1 
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DO  110  1-1 . IL-1 

XL(I)-(XB(I+IM2)+AL/4.)/AL 
CPL(I)— CP(1+IM2) 

1 10  CONTINUE 

XL(IL)«XU(1) 

CPL(IL)-CPU(1) 

C. .DEFINE  PLOTTING  AREA  AND  HEADING 

CALL  PHYSORO  .75,3.0) 

CALL  AREA2D(XLEN, YLEN) 

CALL  SCLPICISCALE2; 

CALL  YAXANG(0 . 0) 

CALL  XTICKS  (5) 

CALL  YTICKS  (5) 

CALL  GRAF(X1 .DX.X2.Y1 ,DY,Y2) 

CALL  THKFRM( . 02) 

CALL  FRAME 

CALI  RESET ( ' THKFRM ’ ) 

C. .PARAMETERS 

V*B  -  2.65 
HB  -  .25 
XBLK  «  XLEN-3.5 
YBLK  -  YLEN+0 . 4 

C  CALL  BLTREC(XBLK , YBLK.KB ,HB ,0.0,0.015) 

C  CALL  BLKEY(IDI) 

C  CALL  BLOFF(IDI) 

XM  -  XBLK  +  0.  10 

YM  -  YBLK  +  0  05 

C  CALL  MESSAG(V(RESSUR£)  D(  1STRI8UTI0N)$\  100.XM.  YM) 

W8  -  1.55 
HB  -  1.05 
XBLK  -  XLEN-1  6 
YBLK  -  YLEN-1 . 1 

CALL  BLTREC(XBLK. YBLK, WB.HB, 00,0  015) 

CALL  BLKEY( ID2) 

CALL  BLOFF(ID2) 

XM  -  XBLK  +  0.  10 

YM  -  YBLK  +  005 
CALL  HEIGHT (0.15) 

CALL  MESSAG('R(E  “)$' . 100.XM, YM) 

yp— Yljt  S 

CALL  REALNO(RE.-I.XR.YM) 

XM  -  XBLK  +  01 

YM  -  YBLK  +0.30 
CALL  HEIGHT(0. 15) 

CALL  MESSAGC(RF  »)$'  .100.XM.YM) 

XR  -  XM  +  0.6 

CALL  REALNO(RF , 3 . XR , YM) 

XM  -  XBLK  +  01 

YM  -  YBLK  +0.55 


CALL  HEIGHT (0.15) 

CALL  MESSAG(’(T  -)$ • . 100, XM, YM) 


XR  »  XM  +  0. 6 

CALL  REALNO( T 1 , 3 , XR , YM) 

XM  -  XBLK  +  01 

YM  -  YBLK  +0.80 
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CAUL  HEIGHT(0 . 15) 

CALL  MESSAG('\A)  «$’ , 100.XM.YM) 

XR  -  XM  +  0.6 

CALL  REALNO(ALPD , 1 , XR , YM) 

CALL  BLON(IOI) 

CALL  BLON( ID2) 

C. .PLOT  THE  CP  DISTRIBUTION 

CALL  THKCRVf . 02) 

CALL  MARKER( 15) 

CALL  CURVE(XU,CPU. IU.0) 

CALL  MARKER(IB) 

CALL  DASH 

CALL  CURVEIXEXP ,CPEXP( 1 , 1  EXP) ,16,0) 
CALL  RESET 1 " DASH ' ) 

CALL  CURVE (XL. CPL , I L. 1 ) 

CALL  RESET ( ' THKCRV ’ ) 

C  CALL  GRID(2,1) 

CALL  ENDPL(0) 

GO  TO  1000 

2000  END IF 

C  . . END  OF  CP  LOOP.  BEGIN  LOADS  LOOP’ 

IF(IOPT.GE.I)  THEN 
C. . INPUT  FROM  LOAOS 


1-1 

2100  READ( 10, EMV2200)  NT1 .NTPL.T1 .ALPD .RF.RE.CL1 .CD1 ,CM1 ,CP 
IFfNTI . LT.NTS)  GO  TO  2100 
IF(NT1 .GT.NTL)  GO  TO  2200 
WRITE(6, 11)  NT1.T1.ALPD.CL1.CD1.CM1 
11  FORMAT ( '  ‘ , I5.5F8.3) 

ANG( I )— ALPD 
T( I )— T1 
CL(T)-CL1 
CD(I)-CD1 
CM(  I )— CM1 
I-I+1 

GO  TO  2100 
2200  CONTINUE 
IPTS— 1-1 

C.. INNER  LOOP  —  SUBPLOTS 

IF(IOPT.EQ.I)  CALL  PHYSOR(3 . 0 . 7 . 75) 

IF( IOPT . EQ. 2)  CALL  OREL (0. .-3 . 50) 

IF( IOPT . EQ. 3)  CALL  OREL(0 . ,-3 . 50) 

C.  SET  UP  ALPHA-NUMERICS  AND  AXES  NAMES 

CALL  NOBRDR 
CALL  COMP LX 
CALL  BASALF (’ STAND ’ ) 

CALL  MIXALF('L/CSTD' ) 

CALL  MX3ALF ( 1 L/CGR ’ , 1H\) 

CALL  BLOWUP (SCALE) 

C . . LABEL  AND  DRAW  AXES 


IF(IOPT . EQ. 1 )  THEN 


o  o  o  o  o 


CALL  YNAMECL(IFT)  C(OEFFICIENT)$\  100) 

Y1— .5 

Y2-2. 

DY-.5 
END  IF 

I F ( I OPT . EQ . 2 )  THEN 

CALL  YNAME( ’D(RAG)  C(OEFFIC IENT)$ • , 1 00) 
Y1— .2 
Y2-1  . 

DY-.2 
END  IF 

I F ( IOPT. EO.3)  THEN 
CALL  XNAME('\A)$\100) 

CALL  YNAME(  ’M(OMENT)  C(OEFFICIENT)$M00) 
Y1— .15 
Y2-.15 
DY-.05 
ENDIF 

CALL  AREA20(XLEN, YLEN) 

CALL  YAXANG(0. 0) 

CALL  XTICKS(5) 

CALL  YTICKS(5) 

CALL  GRAF(X1 ,DX,X2,Y1 .DY.Y2) 

C . . PARAMETERS 

I F ( I PAR . EQ . 1  .AND.  IOPT.EQ.1)  THEN 
W8  -  1.55 
HB  -  0.55 
XBLK  -  0.05 
YBLK  -  YLEN-0.6 

CALL  BLTREC(XBLK. YBLK. WB.HB. 0.0. 0.015) 
CALL  BLKEYf  ID3) 

CALL  BLOFF( ID3) 

XM  -  XBLK  +  0.10 

YM  -  YBLK  +  0.05 
CALL  HEIGHT(0. 15) 

CALL  MESSAGE R(E  -)$’ . 100.XM.YM) 

XR-XM4-.5 

CALL  REALNO(RE.-1 .XR.YM) 

XM  -  XBLK  +  0.1 

YM  -  YBLK  +  0.30 
CALL  HEIGHT(0. 15) 

CALL  MESSAG( 1 (RF  ■)$• . 100.XM.YM) 

XR  -  XM  +  0.6 

CALL  REALNO(RF.3.XR,YM) 

CALL  BLON( ID3) 

ENDIF 

C. .PLOT  THE  LOADS  DISTRIBUTION 
CALL  THKCRV( .017) 

IF(IOPT.EQ.I)  CALL  CURVE(ANG,CL, IPTS.0) 

IF( IOPT. EO. 2)  CALL  CURVE( ANG .CD , IPTS , 0) 
ir(IOPT.EQ.3)  CALL  CURVE(ANG .CM , IPTS  , 0) 

CALL  DASH 

I F( IOPT.EQ.1)  CALL  CURVE( ALPHA , Cl  EXP , 20.0) 
I F( IOPT . EQ. 2)  CALL  CURVE( ALPHA .CDEXP , 20,0) 
IF(  IOPT. EQ.3)  CALL  CURVE(ALPHA .CMEXP, 20.0) 
CALL  RESET ( ’DASH-) 

CALL  RESET ( ’THKCRV') 

CALL  THKCRV ( . 007 ) 


CALL  GR1D0 .1 ) 

CALL  RESET ( 'THKCRV’ ) 

CALL  THKCRV(.0001) 

CALL  GRID(2,2) 

CALL  RESET ( ’ THKCRV  * ) 

CALL  THKFRM{ .023) 

CALL  FRAME 

CALL  RESET ( * THKFRM  * ) 

FINISH  SUBPLOTS 

CALL  ENOGR(0) 

IOPT-IOPT+1 

IF(IOPT.LE.3)  GO  TO  2290 
CALL  ENOPL(0) 

ENDIF 


»ENO  OF  LOADS  LOOP* 


C.  PLOTTING  FINISHED 

3000  CALL  DONEPL 
STOP 
END 


m 
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PROGRAM  CARPET 

PURPOSE:  Creotes  a  single  plot  with  multiple  Cp  vs  X/C  curves  in  o 
carpet  plot  format.  The  Cp  and  X/C  scales  are  marked  on  the 
first  curve  drawn.  The  other  curves  are  drawn  to  the  same 
scale,  but  axes  are  not  drawn  although  the  Cp«0.0  level  is 
shown  with  a  dotted  line. 


EXTERNAL  REFERENCES: 
MODULE  DESCRIPT I 
DIP  Initial i z 
DISSPLA  Graphics 
OPENER  Prompts  f 
PLFREE  Writes  "f 
PLTYTL  Writes  a 


DESCRIPTION 

Initializes  file  in  NASA-Ames  Device  Independent  Plot  format. 
Graphics  software. 

Prompts  for  the  name  of  a  sequential  file  and  opens  i t . (PROGTOOLS) 
Writes  "free  values"  with  headings  on  a  plot  (GRAPHLIB) 

Writes  a  centered  character  string  on  the  plot.  (GRAPHLIB) 


DEVELOPMENT  HISTORY: 

DATE  INITIALS  DESCRIPTION 
Nov. 1986  RCL  Original  design  and  implementation. 

AUTHOR:  Rosalie  C.  Lefkowitz  Sterling  Software,  Palo  Alto,  CA 


IMPLICIT  NONE 
•  Parameter  constants: 

INTEGER 

>  LUNCRT,  LUNDAT,  LUNDIP,  LUNIN.  LUNKBD.  MAXI.  MAXI2 ,  MAXJ , 

>  MXSTR.  MAXTIT,  NFREE 
PARAMETER 

>  (  LUNCRT  -  6. 

>  LUNDAT  -  8, 

>  LUNDIP  -  10, 

>  LUNIN  -  12. 

>  LUNKBD  -  5. 

>  MAXI  -  40, 

>  MAX  1 2  -  40.  I  MAX  1/2  +  1  . 

>  MAXJ  -  20. 

>  MXSTR  -  81 . 

>  NFREE  -  3  ) 


C  *  Variables: 


>  ALPHA(MAXJ) ,  CPTH(MAXI.MAXJ),  FREVAL(NFREE)  .  XTH(MAXI), 

>  X2(3) .  Y2(3) 

DIMENSION  WORK (MAX I ) 

REAL 

>  FRELEN.  HIGH.  HITFRE,  HITITL,  WIDE. 

>  XOFSET ,  XOR.  YOFSET .  YOR,  YINCH,  YFREE,  YTITLE. 

>  XMAX.  XMIN,  XSTEP , 

>  CPMAX ,  CPMIN,  CPSTEP 

CHARACTER 

>  DATFI 1*80 , 

>  FRETXT(3) • ( 10) . 


sw 


>  TITLE* (MXSTR) 


INTEGER 

>  I.  J.  NOGFRE(NFREE) 

DATA 

>  CPMIN.CPSTEP.CPMAX/  0. .-20. .-20./, 

>  DATFIL  /  ’ CARPET . DAT '  /. 

>  FRELEN/10.0/. 

>  HIGH/3.0/, 

>  HITFRE/0 . 09/. 

>  HITITL/0. 14/. 

>  NOCFRE/2,  2.  3/. 

>  WIDE/2.8/. 

>  XMIN , XSTEP , XMAX/0 . , 1 . , 1 . / . 

>  X2/0.  ,0.  .  1  ./. 

>  XOFSET/ . 095/ . 

>  XOR/0.0/. 

>  Y  FREE/0. 5/, 

>  YOFSET/ . 1 6/ 

>  YOR/0.0/. 

>  YTITLE/6 . 5/ 

C  *  Free  text  for  plot  labeling  and  label  for  alpha  value  list: 

WRITE(  FRETXT(I),  1100  )  •AMPLITUDE’ 

WRITE(  FRETXTf 2) .  1100  )  ’ST.  MEAN’ 

WRITE(  FRETXT (3) .  1100  )  'RED.  FREQ.’ 

C  •  Initialize  plotting  device: 

CALL  DIP  (  LUNOIP,  ’CARPET .DIP’ .  100  ) 


C  •  Prompt  for  and  open  data  file: 

CALL  OPENER  (  LUNCRT , 

>  'Enter  data  file  name  (defoult  is  CARPET.DAT): 

>  LUNKBD ,  DATFIL.  LUNDAT,  ’OLD’  ) 


C  •  Read  the  data  file: 

READ  (LUNDAT. 1000)  TITLE 
WRITE(LUNCRT . 1 000)  TITLE 
READ  (LUNDAT.*)  (  XTH( I ) , 1-1 .MAXI  ) 

READ  (LUNDAT,*)  (  FREVAL( I ) , 1-1 .NFREE  ) 

DO  10  J»1 . MAX J 

READ  (LUNDAT,*)  ALPHA(J) 

READ  (LUNDAT.*)  (  CPTH( I . J) , 1-1 ,MAXI  ) 
10  CONTINUE 


*v 
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CALL 

RESET 

CALL1) 

CALL 

SETDEV 

(  0.  0  ) 

CALL 

HWROT 

(  ’MOVIE’ ) 

"V 

CALL 

HWSCAL 

( ’NONE’) 

CALL 

N06RDR 

N 

CALL 

GRACE 

(  0-0  ) 

vi 

CALL 

NOCHEK 

CALL 

PAGE  ( 

11.0.  8.5  ) 

?= 

V*v 

CALL 

BASAL F 

(  ’STANDARD’  ) 

CALL 

MIXALF 

(  ’ L/CGREEK ’ ) 

> 

i 


a 


ttBMiiMtailll«^^ 


1 


“1 


C  »  Start  with  title  and  free  text/free  values: 

CALL  PHYSOR  (  XOR,  YOR  ) 

CALL  AREA2D  (  10.0.  8.0  ) 

CALL  PLTYTL  (  0.0,  FRELEN,  YTITLE,  HITITL,  TITLE  ) 

CALL  PLFREE  (  NFREE.  FRETXT,  FREVAL,  NDGFRE.  H1TFRE,  YFREE, 
>  FRELEN  ) 

CALL  ENDGR(0) 


Then  draw  Cp  curves: 

XOR  -  2.0  -  XOFSET 
YOR  -  1.0-  YOFSET 
00  500  J  -1.MAXJ 

XOR  -  XOR  +  XOFSET 
YOR  ■  YOR  +  YOFSET 
CALL  PHYSOR  (  XOR,  YOR  ) 

CALL  AREA2D  (WIDE,  HIGH) 

CALL  CROSS 

IF  (  J.EQ.1  )  THEN 

First  plot  has  full  labeled  axes: 

CALL  XNAME  (  ’X/C' ,  3  ) 

CALL  YNAME  (  '-Cp',  3  ) 

ELSE 

Later  curves  are  drawn  without  labeled  axes: 

CALL  XNAME  \  0  ) 

CALL  YNAME  (’  0  ) 

END  IF 

CALL  GRAF  (XMIN,  XSTEP,  XMAX,  CPM1N,  CPSTEP,  CPMAX) 
CALL  CURVE  (  XTH(1).  CPTH(I.J).  MAXI2,  0  ) 

IF  (  J.GT.  1)  THEN 

Draw  dotted  pseudo-axes: 

Y2(1)  -  CPTH(I.J) 


Y2(2)  -  0.0 
Y2(3)  -  0.0 
CALL  DOT 

CALL  CURVE  (  X2,  Y2 .  3.  0  ) 

CALL  RESET  (  'DOT1  ) 

END  IF 

CALL  HEIGHT  (  HITFRE  ) 

C  •  How  many  inches  from  the  origin  is  Y-  0  ? 

YINCH  -  YPOSN  (  0.0,  0.0  ) 

C  •  Now  print  the  value  of  ALPHA: 

CALL  REALNO  (  ALPHA(J).  2,  WIDE  +0.2,  YINCH  ) 

C  •  Identify  the  alpha  value  list: 

IF  ( J . EQ.MAXJ  ) 

>  CALL  MESSAG  (’(a)  deg.$’,  100.  WIDE  +  0.2.  YINCH  +  0.2  ) 
CALL  ENDGR  (0) 

500  CONTINUE 

C  •  Termination: 

CALL  ENOPL  (  0  ) 

CALL  DONEPL 

STOP  '  Normal  termination.  Plot  file  is  in  CARPET. DIP. 

1000  FORMAT  (A) 

1100  FORMAT  (A10) 


DISCLAIMER 


The  following  is  o  site-specific,  step-by-step  example  intended  for  the 
novice  user  to  employ  the  Wu  code  at  the  NASA-AMES  computational  focility. 

It  makes  no  attempt  to  be  all  encompassing,  but  rather  tries  to  provide  an 
adequate  amount  of  information  to  shorten  the  learning  period  required  to 
gain  a  minimum  working  knowledge  of  the  facilities.  Much  of  it  can  be  applied 
to  the  use  of  other  codes  at  AMES,  but  it  will  have  only  limited  application 
at  other  sites. 

There  are  a  number  of  useful  publications  which  can  be  helpful  in 
utilizing  the  NASA-AMES  computational  facilities.  Included  are: 

t.  INTRODUCTION  TO  VAX/VMS  AT  AMES 

2.  CrayVAX  DECnet  Interface  User's  Guide 

3.  A  GUIDE  TO  GENERATING  MOVIES  USING  PL0T3D  AND  GAS 

4.  PL0T3D  and  PL0T3X  Version  3.5  (3D  for  VAX  and  3X  for  Cray) 

Joan  Thompson  and  Roseal ie  Lefkowitz  usually  have  these  publications. 

Their  office  is  in  Bldg  227,  Rm  102.  Both  ladies  are  very  helpful  and 
generous  with  their  time. 

You  will  need  access  to  a  VAX  account  and  required  password.  See  an 
NPS  representat i ve  for  this  information.  'Jianps'  is  on  occount  on  RALph 
currently  available  for  NPS  use,  which  is  convenient  as  DISSPLA  and  PLOT3D 
are  installed  on  RALph.  Individual  VAX  are  also  referred  to  as  'nodes’. 

The  following  description  assumes  the'user  is  at  a  graphics-capable 
terminal  (or  'green  screen').  Use  for  standard  terminals  will  be  identical 
with  the  exception  of  plotting  locations  available.  Remote  use  does  not 
allow  for  full  screen  editing. 

First,  logon  by  typing  'c  •••'  ,  where  'c’is  for  connect,  and  is  for 

the  first  three  letters  of  the  node  you  are  logging  onto.  You  will  then  be 
prompted  for  your  password  (pw).  Upper  or  lower  case  is  allowable  for  almost 
a  I  I  commands. 

All  regular  keyboard  entries  are  submitted  to  the  VAX  via  the  carriage 
return  <cr>.  with  the  exception  of  command  line  entries  which  submitted  via 
the  'Enter'  key. 

Now  that  you  ore  logged  into  the  system,  you  can  check  to  see  what  you  have 
immediate  access  to  by  typing  'dir'. 

If  the  file  reeded  is  in  a  subdirectory,  type  ' sd  [•••]',  where  '•••'  is 
the  subdirectory  name. 

Type  'dir'  (or  'list'  for  added  file  information)  again. 

To  edit  or  to  look  at  the  file,  type  'edt  fn’  .where  fn  is  the  filename, 
(i.e.  this.fi  I e ; 1  or  my.dat)  The  default  version  (or  number  after  the 
semicolon)  is  the  highest  number. 

You  are  now  in  the  line  editing  mode.  To  go  to  full  screen  or  keypad 
editing,  type'c'  and  <cr>. 

The  keypad  is  the  at  the  right  end  of  the  keyboard.  See  the  "Introduction  to 
VAX/VMS  at  AMES”  for  more  information.  AMES  will  eventually  be  converting  to 
UNIX  (which,  like  “VAX"  is  a  unique  operating  system)  and  have  some  different 
functions  and  options  than  are  on  "DEC".  DEC  and  VAX  ore  normally  used 
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i nt e rchongab I y .  but  OEC  is  propsrly  the  name  of  the  manuf ac t u re r , wh i I e  VAX 
is  the  operating  system  itself 

After  you  are  finished  changing  or  looking  at  the  file,  enter 
'PF1'  and  then  keypad  '7'.  This  w  I l  put  you  in  the  command  line  mode. 

To  save  any  changes  and  create  a  new  version  of  the  file,  type  'exit'  enter 
If  no  changes  have  been  made  or  you  don't  want  to  save  the  changes  made,  type 
'quit'  and  ' Enter ' 

To  submit  a  job  to  the  CRAY  X/MP-A8,  type  'csub'  and  the  fn.  You  will  then 
be  prompted  for  your  pw.  For  additional  information,  see  "Cray  VAX  DECnet 
Interface  User's  Guide". 


The  following  is  an  example  of  a  series  of  Cray  and  VAX  runs  to  employ  the 
Wu  code  and  plotting  options. 


NOTES  XI LL  BE  INSERTED  THROUGHOUT. 


FOR  CRAY  JOB  CONTROL  LANGUAGE  (JCL).AN  '*.'  KILL  COMMENT  OUT  A  LINE. 

ALL  CRAY  JCL  LINES  MUST  END  WITH  A  PERIOO. 

THE  FIRST  LINE  OF  THE  JCL  MUST  BE  THE  ' JOB— '  LINE.  DON'T  HAVE  A  BLANK 
LINE  FOR  THE  FIRST  LINE  OR  THE  FIRST  SPACE  ON  THE  FIRST  LINE!  IT  WON'T 
WORK. 

INSERT  YOUR  OWN  JOB  NAME  (JN),  USER  ID  (US)  ANO  USER  PASSWORD  (UPW) . 

YOUR  JN  MAY  BE  REPEATED  FOR  THE  NAME  AND  VAX  ADDRESS  (VAX  OR  NODE 
AODRESS  FOR  THIS  CASE  IS  FML.  OTHERS  INCLUDE  RALph.MARs.TOM.  etc.) 

•X-'  IS  FOR  YOUR  AMES  PHONE  EXTENSION. 

AS  OF  19APRIL88,  THE  FIRST  TWO  TIME  LIMIT  SIZES  ARE  T-150  (CPU  SECONDS) 
AND  T-1800.  THE  FIRST  WILL  BE  ADEQUATE  FOR  ALL  JOBS  UP  TO  ABOUT 
100  TIME  STEPS.  THE  CPU  LENGTH  OF  EACH  TIME  STEP  tS  DIRECTLY  PRORTIONAL 
THE  SIZE  OF  THE  INVISCID  REGION  BEING  CALCULATED.  THEREFORE,  PRE-STALL 
ANGLE  OF  ATTACK  CALCULATIONS  WILL  BE  CLOSER  TO  .5  CPU  SECONDS  PER  TIME 
STEP  AND  FULLY  STALLED  CONDITIONS  WILL  REQUIRE  1  TO  1 . 3  CPU  SECONDS 
PER  TIME  STEP. 


THE  FOLLOWING  IS  A  POSSIBLE  SERIES  OF  COMPUTER  RUNS. 


RUN  1 :  SAVE  DATA  FOR  STEADY  STATE  CASE 


JOB,  JN-*****  .T-150. 

ACCOUNT , AC—* • • * • . US—* *  * • • , UPW-* 


X— 4269  FML: : < 


ACCESS. Dt**SENOVAX ,  ID-STTRDM.OWN-RFTRDM. 

•  . 

*.  Compile,  load  and  run 

•  . 

*.CFT.ON-CSTAX.(USE  THESE  OPTIONS  FOR  ADDITIONAL  DEBUGGING  INFORMATION) 
CFT.ON-A.OFF-CST. (USE  THESE  OPTIONS  FOR  ALL  REGULAR  RUNS.  THEY  PROVIDE  A 
LDR.  MINIMUM  OF  EXTRA  OUTPUT  INFORMATION.) 

•  . 

NOTE:  ENTER  YOUR  DESIRED  DIRECTORY  FOR  THE  OUTPUT  DATA  TO  BE  SENT  TO. 

SENDING  THEM  TO  SCRATCH  WILL  HELP  AVOID  EXCEEDING  YOUR  ALOTTED 
DISC  QUOTA  SPACE.  YOUR  JN  CAN  BE  USED  FOR  THE  •*•  BELOW. 

SENOVAX.DN-FT01 .VON-' DUA1 :[ SCRATCH. ••«] TAP El .OAT’ .  (UNCOUJENT  THESE 
SENDVAX ,DN— FT03 , VDN-' DUA1 :  SCRATCH .*••  TAPE3 . DAT ' .  LINES  WHEN  USING 
SENOVAX. DN-FT 1 4.  VDf*. 'DUA1  .  [SCRATCH.  •*•  JTAPE14.DAT  ’  .  THE  DISSPLA  ROUTINES) 

NOTES:  DISSPLA  ROUTINES  INCLUDE  PLOT1.PLOT2  AND  PLOT3. 

UNCONMENT  THE  FOLLOWING  LINE  WHEN  GENERATING  DATA  FOR  PLOT3D . 

SENDVAX.DN— FT10,  VDfW*'  FAR0 :  [  .  NAMEJGRID.DAT'  . 

*  , 

REWIND .DN-FT02 . 

SK I PF.Df^S IN. (BYPASSES  ANY  SPURIOUS  LINES  AFTER  DATA.) 

*  . 

NOTE:  THE  PON  IS  FOR  THE  DATA  SAVED  FOR  EITHER  THE  INITIAL  STEADY  STATE  CASE 
(WHICH  IS  NEEDED  FOR  ALL  DYNAMIC  CASES)  OR  FOR  RESTARTS. 

ACCESS. DN-FT07. PON-TEST  1 . 

•  . 

NOTE:  THIS  IS  THE  START  OF  THE  JCL  FOR  THE  SECOND  JOB  TO  BE  RUN.  SUBMITTING 
MULTIPLE  JOBS  IS  IN  THIS  MANNER  IS  CALLED  'JOB  CHAINING.' 
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» .CFT.ON-CSTAX. 

CFT.ON-A.OFF-CST. 

LDR . 

•  . 

SAVE, OhH>FT08. PDN-TEST1  .  (USE  ONLY  WHEN  THE  DATA  GENERATED  WILL) 

.  .  (BE  NEEDED  FOR  A  SUBSEQUENT  RUN.) 

S ENDVAX  , DN«FT 09  ,  VDN" ’ DUA 1  : [SCRATCH. 1TAPE9 , DAT  *  . 

SENDVAX , DN-FT04 , VDN“ ’ DUA 1 : [ SCRATCH . • • • j  T APE4 . DAT ‘ . 

SENOVAX.ON-FT1 1 , VDN“ ‘ [ . JQ. DAT ' . 

/EOF 

0  .15  5.0  1  (ICST-0  TO  COMPUTE  THE  STEADY  STATE  CASE.) 

.0001  .0005  .0005  75  100 
.4  .3  .6 

10  35  50  74  45  40 
.10  .10  25 

25  25  25  (USE  THESE  VALUES  TO  ALLOW  A  CHECK  ON  THE  OUTPUT  SOLUTION.) 


ICST , RF.ALPS , ICTUR 

WM I N . DFVX , DRMX . KMAX , NCC 

URBI ,URBP,URR 

IV1 . IB1 , IB2, IV2.NPL.NLB 

DTI , DTINC.NTMAX 

NTPL,  NTOUT ,  NTL.O 


RUN  2:  A  SMALL  NUMBER  OF  TIME  STEPS  TO  CONFIRM  THE  RESULTS  ARE  AS  DESIRED. 
NOTE:  ALL  THE  ABOVE  JCL  IS  THE  SAME  EXCEPT: 

*. SAVE. DhN.FT08. PDN-TEST1  .  (USE  ONLY  WHEN  THE  DATA  GENERATED  WILL) 

.  (BE  NEEDED  FOR  A  SUBSEQUENT  RUN  ) 

AND 

•  .  SENDVAX  ,DtA*FT1 0  ,VD(^‘  FAR0  :  [  .  NAME]GRID  ,  DAT  ‘  . 

4  .15  5.0  1  (ICST-1,2.3  OR  4  AS  DESIRED.) 

.0001  .0005  .0005  75  100 
4  3  6 

10  35  ,0  74  45  40 
.10  10  25 

25  25  25  (USE  THESE  VALUES  TO  ALLOW  A  CHECK  ON  THE  OUTPUT  SOLUTION.) 


ICST. RF.ALPS. ICTUR 

WM I N , DFMX , DRMX , KMAX . NCC 

URBI .URBP.URR 

IV1 . IB1 . IB2 . IV2 . NPL.NLB 

DTI .DTINC.NTMAX 

NTPL. NTOUT, NT LO 


RUN  3:  THF  FINAL  DATA  CAN  BE  SAVED  ON  TAPE  IF  A  RESTART  IS  TO  BE  USED. 
NOTE:  ALL  THE  ABOVE  JCL  IS  THE  SAME  EXCEPT 
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SAVE. DN-FT08, P0N-TEST1 . (USE  ONLY  WHEN  THE  DATA  GENERATED  WILL) 

. .  (BE  NEEDED  FOR  A  SUBSEQUENT  RUN.) 

A  .15  5.0  1  (ICST-1.2.3  OR  4  AS  DESIRED.) 

.0001  .0005  .0005  75  100 
.4  .3  .6 

10  35  50  74  45  40 
.10.10  500 

50  50  50  (SELECT  OUTPUT  TIME  STEP  VALUES  AS  DESIRED.  THESE  VALUES 

WILL  GIVE  OUTPUT  AND  PLOTTING  DATA  IN  THE  SAME  4.0  SECOND 
INTERVALS  AS  THE  MANUAL  DOES.  ALTERNATIVELY,  ENTER  THE 
SPECIFIC  ALLP  (OR  ALPHA  IN  DEGREES)  VALUES  IN  ZONST. 

C OAK  ENT  OR  UNCOMMENT  THE  LINES  ASSOCIATED  WITH  ALL P 
DEPENDING  ON  THE  MANNER  OF  CHOOSING  THE  OUTPUT  CONTROL.) 


ICST , RF . ALPS , ICTUR 

WM1N.DFMX.DRMX.KMAX.NCC 

URB1 .URBP.URR 

IV1 , IB1 . IB2, IV2.NPL.NLB 

DTI , DTINC.NTMAX 

NTPL . NTOUT.NTLO 


AFTER  EACH  RUN.  LOOK  AT  THE  OUTPUT  FILE  TO  SEE  IF  THE  JOB  HAS  RUN  TO 
COMPLETION.  THIS  IS  INDICATED  IF  THE  CPU  TIME  IS  REASONABLE  FOR  THE  NUMBER 
OF  TIME  STEPS  RUN.  IF  THERE  ARE  NO  MESSAGES  THAT  THE  JOB  WAS  ABORTED.  AND 
THAT  NORMAL  DATA  TRANSFER  WAS  ACCOMPLISHED.  THEN.  IF  IT  HAS,  PLOTTING  CAN 
BE  DONE. 


IF  THE  WU  PLOTTING  ROUTINES  (OISSPLA)  ARE  TO  BE  USED  ENSURE.  THE  FOLLOWING 
ARE  AVAILABLE  IN  YOUR  DIRECTORY. 

FOR  PLOT1 : 

1 .  PLOT  1 .EXE 

2.  PLOT 1  DAT 

3.  PLOT1.COM 

IF  YOU  DON'T  HAVE  A  VERSION  OF  PL0T1.EXE  IN  YOUR  DIRECTORY.  COMPLETE 
THE  NEXT  TWO  STEPS. 

COMPILE  PLOT 1 . FOR  BY  TYPING: 

FOR  PL0T1  (FOR  AND  HIGHEST  VERSION 

NUMBER  ARE  THE  DEFAULT.)  THIS  WILL  CREATE  PLOT  1. OBJ 

LINK  PL0T1  OBJ  BY  TYPING: 

LINK  PL0T1 ,SYS$LIBRARY: INTL1B/LIB , DI SSPLA/LIB, I NT/LIB 
THIS  WILL  CREATE  PL01  EXE  THIS  IS  THE  FILE  NEEDED  TO 
ACTUALLY  DO  THE  PLOTTING. 

EDT  PLOT1.COM  TO  ENSURE  THAT  THE  PROPER  NAMES  ARE  INCLUDED  IN  THE  FILE  NAMES. 
THESE  SHOULD  BE  THE  SAME  NAMES  AS  IN  YOUR  CRAY  JCL. 

YOU  SHOULD  ALSO  HAVE  SAVED  THE  DATA  BY  UNCOMMENT  I NO  THE  APPROPRIATE  LINES 
IN  YOUR  JCL.  SEE  ABOVE. 

PLOT 1 .COM  CAN  BE: 

$  GRAPHICS 

$  IF  " ' ' F$MOOE( ) ' "  EQS.  "BATCH"  THEN  SET  DEFAULT  FARO : [ J IANPS . •• ] 

$  DEFINE/USER  FOR001  DUA1 : [SCRATCH .•«. JTAPE1 . DAT 
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$  DEFINE/USER  FOR005  PLOT  1  . l)AT 
$  RUN  PLOT1 

WHERE  •  IS  YOUR  JN  AND  ••  IS  YOUR  DIRECTORY 

TO  CREATE  THE  PLOTS  AVAILABLE  FROM  PLOT1 ,  TYPE: 
OP LOTI 


FOR  PLOT2: 


1.  PLOT2.EXE 

2.  PLOT2.DAT 

3.  PL2.COM 


IF  YOU  DON'T  HAVE  PL0T2.EXE  FOLLOW  THE  ABOVE  DIRECTIONS  FOR  PL0T1 


PLOT2.COM  CAN  BE: 


GRAPHICS 

IF  " ’ ’ F$MOOE( ) ' "  . EQS 
DEFINE/USER  FOR001  DUA1 
DEFINE/USER  FOR003  DUA1 
DEFINE/USER  FOR004  DUA1 
DEFINE/USER  FOR009  DUA1 
DEFINE/USER  FOR014  DUA1 
DEFINE/USER  FOR005  PL0T2.DAT 
RUN  PL0T2 

DIPQMS/DELETE  PL0T2.DIP 


BATCH"  THEN  SET  DEFAULT  FAR0: [JIANPS. *•] 


SCRATCH. 
SCRATCH. 
SCRATCH. ••• 
SCRATCH. ••• 
SCRATCH. ••• 


TAPE! . DAT 
TAPE3.DAT 
TAPE4 . DAT 
TAPE9 . DAT 
TAPE14.DAT 


TO  CREATE  THE  PLOTS  AVAILABLE  FROM  PL0T2 ,  TYPE: 
•PL2 


FOR  PLOTS: 

1.  LOADS . DAT 

2.  LOADS.EXE 

3.  LOADS.COM 

4.  PL0T3.EXE 

5.  PL0T3.DAT 

IF  YOU  DON'T  HAVE  PLOTS- EXE  OR  LOADS.EXE.  SEE  ABOVE. 

LOADS.COM  CAN  BE: 

$  GRAPHICS 

$  IF  ” ’ ’ F$MODE( ) ' "  .EQS.  "BATCH"  THEN  SET  DEFAULT  FARO: [JIANPS. •* . • 
%  DEFINE/USER  FOR003  DUA1 : [SCRATCH . •♦•lTAPES . DAT 
$  DEFINE/USER  FOR004  DUA1 : [SCRATCH . •••JTAPE4 . DAT 
$  DEFINE/USER  FOR005  LOADS . DAT 
$  RUN  LOADS 

$  !  OEFINE/USER  FOR024  NACA.DAT 
$  DEFINE/USER  FOR010  FORO10.DAT 
$  DEFINE/USER  FOR014  DUA1 : [SCRATCH . •••]TAPE14 . DAT 
$  DEFINE/USER  FOR005  PLOTS. OAT 
$  RUN  PLOTS 


TO  CREATE  THE  PLOTS  AVAILABLE  FROM  PLOTS,  TYPE: 
OLOAOS 


FOR  USING  PL0T3D ,  THESE  ARE  A  SET  OF  POSSIBLE  DATA  ENTRIES: 


GRAPHICS 

PLOT3D 

READ/20 

THE  FOLLOWING  ARE  RESPONSES  TO  PROMPTS: 

GRID 

Q 

SUBSETS 
1  80  5 

1  60  5 

1 


MINMAX  -55-55 
FU  200 
WALLS 
1  80 


1 


1 


VECTORS 

Y 

.1 

LINES 

Y 
.5 


P/2D/TEK 

THIS  WILL  GIVE  A  REPRESENTATIVE  PLOT  OF  THE  VELOCITY  VECTORS.  BY  CHANGING 
THE  MINMAX  AND  SUBSETS  VALUES,  DIFFERENT  ASPECTS  OF  THE  FLOW  CAN  BE 
HIGHLIGHTED. 

THE  LAST  COMMAND  WILL  GENERATE  A  PLOT  ON  THE  TERMINAL.  IF  A  HARD  COPY 
IS  DESIRED,  AFTER  THE  PLOT  IS  VIEWED  OR  IF  YOU  ARE  SURE  YOU  KNOW  WHAT  THE  PLOT 
WILL  LOOK  LIKE.  TYPE:  9/ DIP.  DO  THIS  FOR  EACH  OF  THE  PLOTS  THAT  YOU  WISH 
TO  SAVE.  WHEN  FINISHED  WITH  PLOT 3D,  EXIT  FROM  IT  AND  TYPE:  DIPQMS  Q.  THIS 
WILL  QUEUE  YOUR  PLOTS  TO  THE  LASER  PRINTER  ON  RALPH.  USE  THE  SAME  TECHNIQUE 
TO  PRODUCE  HARD  COPIES  OF  THE  DISPLAY  PLOTS.  THERE  IS  ALSO  A  BATCH  SUBMITTAL 
FILE  FOR  PLOT2  AND  PLOT3 . 
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